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High Fidelity CFD Analysis and Validation of Rotorcraft Gearbox 
Aerodynamics Under Operational and Oil-Out Conditions 

Final Report 


Robert F. Kunz 

The Pennsylvania State University 
University Park, Pennsylvania 16804 


Executive Summary 

This report, submitted by the Pennsylvania State University Applied Research Laboratory, 
represents a final deliverable for NASA Cooperative Agreement NNX07AB34A, “High Fidelity CFD 
Analysis and Validation of Rotorcraft Gearbox Aerodynamics under Operational and Oil-Out 
Conditions.” 

This report is submitted with the two other final deliverables of the contract: 1) NPHASE-PSU 
Theory, User’s and Reference Manual Version 3.16, and, 2) NPHASE-PSU V3.16 source code, 
installation scripts, tutorial and instructions transferred electronically (CD contains this report, the 
manual and nphase-05-31-2013.tar.gz) 

This report summarizes the many technical accomplishments achieved under the two-part NRA. 
The technical and other accomplishments of the NRA are summarized here: 

1. The CFD code, NPHASE-PSU, was instrumented to accommodate the numerous challenges 
involved in modeling gearbox systems, including: 

multiphase flow modeling for oil droplets, oil films and compressible air 
gear motion in the relative and absolute frames of reference 
gear tooth contact (in the context of overset meshing) 
conjugate heat transfer 

2. These NPHASE-PSU capabilities were brought to bear in analyzing gearbox aerodynamics 
in the context of other modern CFD practices including high order discretization, parallel 
programming/execution (on NASA NAS HPC resources), implicit multi-component 
numerics, and 2-4 equation/LES turbulence modeling. 

3. Many verification and validation cases were carried out 

4. Significant new understanding of the physics of gearbox windage was achieved. This has led 
directly to design guidance (including a patent and several rotorcraft industry sponsored 
projects) 

5. This research has led to five conference publications 1 ’'', two journal publications vl ’ vn , two PhD 
students/theses vllwx , and one (provisional) patent x . 

i) Hill, M.J., Kunz, R.F., Medvitz, R.B., Noack, R.W., Long, L.N., Morris, P.J., Handschuh, 

R.F., ’’Application and Validation of Unstructured Overset CFD Technology for 
Rotorcraft Gearbox Windage Aerodynamics Simulations,” American Helicopter 
Society 64th Annual Forum, Montreal, Canada, April, 2008. 

ii) Hill, M.J., Kunz, R.F., Noack, R.W., Handschuh, R.F., Long, L.N., Noack, R.W., ”CFD 

Analysis of Gear Windage Losses: Validation and Parametric Aerodynamic Studies,” 
American Helicopter Society 66th Annual Forum, Phoenix, AZ, April, 2010. 

iii) Kunz, R.F., Hill, M.J., Schmehl, K.J., McIntyre, S.M., “Computational Studies of the 

Roles of Shrouds and Multiphase Flow in High Speed Gear Windage Loss,” 

American Helicopter Society 68th Annual Forum, FortWorth, TX, April, 2012. 
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iv) McIntyre, S.M., Yu, Q., Saribay, Z.B., Kunz, R.F., ’’Fluid and Thermal Dynamics 

Simulation of Spur Gears Under Steady State Operation,” Manuscript accepted for 
American Flelicopter Society 69th Annual Forum, Phoenix, AZ, April, 2013. 

v) Noack, R. W., Boger, D. A., Kunz, R. F., and Carrica, P. M., "Suggar++: An Improved 

General Overset Grid Assembly Capability," AIAA Paper 2009-3992, 19th AIAA 
Computational Fluid Dynamics Conference, San Antonio, TX, June 2009. 

vi) Hill, M.J., Kunz, R.F., Medvitz, R.B., Handschuh, R.F., Long, L.N., Noack, R.W., 

Morris, P.J. “CFD Analysis of Gear Windage Losses: Validation and Parametric 
Aerodynamic Studies,” ASME Journal of Fluids Engineering , 133:03 11 03-1-10, 
March 2011. 

vii) Hill, M.J., Kunz, R.F., Noack, R.W., Long, L.N., Morris, P.J., Handschuh, R.F., ”CFD 

Technology for Rotorcraft Gearbox Windage Aerodynamics Simulations,” Gear 
Technology, 48-55, August, 2009. 

viii) Hill, M.J., A Computational Investigation of Gear Windage , PhD Thesis, Department 
of Aerospace Engineering, The Pennsylvania State University, 2010. 

ix) McIntyre, S.M., Multiphase Modeling of Droplet-Film Flows , PhD Thesis, Department 

of Aerospace Engineering, The Pennsylvania State University, anticipated 2015. 

x) Kunz, R.F., Medvitz, R.B., Hill, M.J., “High Speed Gear Sized and Configured to Reduce 

Windage Loss,” Provisional Patent Application No. 61/477,666. 

6. The NRA project has led to MANY synergistic and follow-on activities, which is a testament 
to the value of this effort: 

Three Bell sponsored VLC projects (helicals) 

Boeing sponsored IRAD effort (spiral bevel) 

VLRCOE LOL project 

Pratt and Whitney Gearbox Lubrication Facility 
Current NASA Glenn instrumentation opportunity 
Benefit of multiphase modeling to Icing NRA with Bell 

7. The PhD student who was funded by the NRA graduated (Matthew Hill) in 201 1 and now 
works at Bell Helicopter 
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Final Report 


In this project, CFD modeling was evolved, applied and validated for gear windage 
aerodynamics. This final report presents the accomplishments carried out under this NRA project. The 
author and his students and colleagues executed research under the NRA from 2007 through 2012, and 
this led to an AHS Forum paper in 2008, which was solicited for and appeared in Gear Technology, an 
AFIS Forum paper in 2010, which was submitted and published in the ASME Journal of Fluids 
Engineering, and an AHS Forum paper in 2012. Additionally some of our research developed under the 
NRA appeared in a paper presented at the 2009 AIAA CFD Conference. 

These four papers represent all of the research carried out under the NRA. We proceed 
chronologically here. 

In the first two years of the program, we enabled NPHASE-PSU with the instrumentation 
required to predict incompressible, single phase (air-only), isolated gear windage simulations, using 
overset meshes, solving on very large meshes in the absolute frame of reference. We demonstrated a 
very convergent dual-time moving mesh numerical scheme, and excellent parallel scaling for the code 
on the large NASA clusters used to perfonn the simulations. The code was verified and validated for 
several canonical flows, and first-of-their-kind fully 3D windage validation cases were carried out 
showing good agreement between CFD and experiment. Viscous/turbulence modeling was identified as 
a shortcoming since loss power was consistently underpredicted for viscous-drag-only spinning disk 
cases. Low Reynolds number (sublayer resolved) turbulence modeling exhibited improved perfonnance 
for these cases. High Reynolds number two equation modeling proved appropriate for modeling the 
moderate speed Diab spur gear suite, due apparently to the dominance of the pressure torques on spin- 
down. The budgets of viscous and pressure components of spin-down torque suggested that viscous 
effects will become much more important, perhaps exceeding 50 percent, for gears with pitchline speeds 
approaching twice that investigated by Diab and the present authors. Some interesting secondary flow 
and pressure gradient physics were identified, although we did not interpret the physics in this early 
work enough to draw any significant conclusions. The 2008 AHS Forum paper which summarizes this 
early research is included here in its entirety. 
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An unstructured overset moving mesh CFD method is adapted, validated and applied to spinning gear 
systems with emphasis on predicting windage losses. Several spur gears and a disc, spinning in air at rotation 
rates up to 1200 s -1 are studied. It is observed that the CFD simulations return good agreement with measured 
windage power loss, as determined by the deceleration of the gears due to torques exerted by viscous and pressure 
forces on the gear surfaces. Hirbulence modeling choices, the relative importance ofviscous and pressure torques 
with gear speed, and the physics of the complex 3D unsteady flow field in the vicinity of the gear teeth are studied. 
The capabilities and challenges associated with the overset mesh approach for enclosed gear train applications 
are demonstrated. 


Introduction 

Gearbox windage refers to the power losses associated 
with rotational deceleration torques exerted on spinning 
gears by aerodynamic forces (pressure and viscous) within 
the air-oil atmosphere present within a gearbox. Windage 
losses are a source of significant heating and fuel consump- 
tion in rotorcraft and other VTOL systems. The weight and 
packaging constraints inherent in these systems require the 
gearing components to be lightweight and heavily loaded. 
The components are also required to operate at high ro- 
tational speeds where windage losses become significant 
with respect to other gearbox losses (rolling, sliding, and 
lubrication losses.) Windage losses are relevant to aircraft 
design for several reasons: 1) They can consume several 
percent of the transmitted power. This has significant im- 
plications for onboard oil cooling requirements and lube 
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system capacity thereby compromising range and standby 
military readiness for rotorcraft and carrier-based aircraft. 
2) Rotorcraft platform survivability under transmission oil- 
out conditions is exacerbated by windage losses which are 
manifested as added dissipative flow heating to these al- 
ready critically thermally stressed systems. 

Despite this significant relevance, design efforts 
throughout the gearing industry aimed at reducing windage 
losses have generally fallen into the trial and ernor category. 
Nevertheless, it has been shown by Winfree 1 and others 
that modest geometric modifications to control the air flow 
path, such as shrouding and baffle configurations, can very 
significantly reduce both windage losses and lubricating 
oil consumption (80% and 40% reductions observed, re- 
spectively, in Ref. 1). However, these hardware specific 
approaches are empirical, expensive, and time consuming, 
and to be relevant need to be performed late in the design 
cycle. 

A host of experimental studies have appeared in the lit- 
erature. 1-8 These studies employ either closed loop sys- 
tems 2-6 or treat isolated gears 1,7,8 where windage (and 
other) losses are determined by measuring spin-down ve- 
locities once the gear+shaft assembly is disconnected from 
the drive torque. Figure 1 shows a diagram of the high 
speed helical gear train test facility at the NASA Glenn 
Research Center. It is a closed-loop system that has been 

5 


setup to study the thermal behavior of aerospace-quality 
gear components under various speeds, loads, and lubri- 
cant flow rates. 5 Generally, these studies parametrize gear 
geometry elements, rotational speed, enclosure geometry, 
and lubrication system characteristics (flow rate, jet loca- 
tion, lubricant rheology), and use dimensional analysis to 
develop correlations for the power losses. These correla- 
tions, although useful in the design process, are inherently 
limited by the large number of system variables and the 
attendant limited range of their - applicability. Of particular 
concern here is the paucity of data/correlations available for 
high speed gears of interest. Indeed, compressibility effects 
are mostly not even considered in the literature, although 
high speed gears can have tip Mach numbers reaching 0.75. 


Slave gearbox ^ -Rotating torque actuator 



Fig. 1 NASA High-Speed Helical Gear Train Facility. 5 

Unfortunately, the physics of these systems are so com- 
plex that to date there have been no attempts made to em- 
ploy many modem elements of 3-D computational fluid 
dynamics (CFD) in analyzing gearbox windage. Recent 
2-D studies 9 were performed using the commercial CFD 
solver FLUENT, where a side correlation factor was used 
to to account for 3D effects (although these authors state 
that work is under way to extend their simulations to 3-D 10 ) 
Specifically, the fluid mechanics involve complex separated 
air flow, disperse multiphase flow (oil droplets) and contin- 
uous multiphase films (lubricating oil on gears), moving 
boundaries in contact, and all modes of heat transfer. Ac- 
cordingly, for a CFD tool to resolve all of the relevant 
physics of this problem it must: 

• support moving meshes (either adaptive/deforming or 
overset) necessary to resolve the gears in relative mo- 
tion and contact at the gear face 

• contain non-equilibrium multiphase flow capability 
(separate continuity and momentum equations for 
each phase; slip between phases) to accommodate the 
disperse mist/droplet and continuous film flows 

• support suitable turbulence modeling to accommo- 
date the complex high speed separated flow within the 
gearbox, and to accurately represent the cascade of en- 
ergy through turbulence scales into viscous heating 


• posses suitable preconditioned time-stepping algorith- 
mics to efficiently accommodate Mach numbers rang- 
ing from near zero through high subsonic 

• support conduction (gears themselves meshed inter- 
nally - solve conduction problem in these elements), 
convection and (for near-failure conditions) radiation 
modeling 

• be parallelized to run efficiently on modem High Per- 
formance Computing (HPC) systems. 

The overall goal of this research is to adapt, validate and 
apply a CFD tool that supports these basic gear - windage 
modeling requirements. This paper represents the first pub- 
lication of the authors’ efforts in this area, and focuses 
on the single phase aspect of the problem, although, the 
CFD code used supports the multiphase flow requirements 
that will be needed in the future. The paper is organized 
as follows: First, the governing equations and numerical 
procedures, including overset grid technology are reported. 
Second, several verification cases are reported, where re- 
sults are compared to analytical solutions for rotating Cou- 
ette flow, flow near an infinite spinning disc, and thermal 
Couette flow. Next several validation cases are reported, 
including windage loss predictions compared to data ob- 
tained by Diab et al., 8 for a spinning disk and four spur 
gears. Gridding requirements and turbulence modeling are 
discussed. 

Theoretical Formulation 

Governing Equations 

The conservation of mass, momentum and energy can be 
written in integral conservation law form for a compress- 
ible flow through a moving mesh as: 

Y t JpdV+Jp(V-W)-dS = 0 (1) 

v s 

p VcN+ I pV(V-W)-dS= - J pdS + j x-dS 
vs s s 

( 2 ) 

j p EdW+ j p H(V-W)-dS= J (Vv) dS+W f +q H 

vs s 

( 3 ) 

In equations 1-3, V is the velocity vector, and W is the 
velocity of the surface element, dS, both in the absolute 
frame of reference. In the present work, all verification and 
validation flows studied are either incompressible, or have 
maximum local absolute Mach numbers of less than 0.35. 
Accordingly, for all simulations presented in this paper, 
an incompressible assumption is invoked and the energy 
equation is not solved (except for the thermal Couette flow 
simulation where it is solved subject to a constant density 
constraint). 
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A high-Reynolds number k-e turbulence model and a 
sublayer resolved hybrid k-e/k-© turbulence model due to 
Menter 11 are used in the studies that follow. No explicit 
transition model was employed as justified, for now, by the 
small contribution of near- axis viscous torques on windage 
loss. 

CFD Numerics and Code 

The CFD code used in this work, NPHASE-PSU 12 
is a parallel face-based, cell-centered, arbitrary- element 
unstructured multiphase flow solver which has been in- 
strumented with overset mesh capability. The baseline 
algorithm follows established segregated pressure based 
methodology. A collocated variable arrangement is used 
and a lagged coefficient linearization is applied. 13 Di- 
agonal dominance preserving, finite volume spatial dis- 
cretization schemes are used for the scalar transport equa- 
tions. Continuity is introduced through a pressure cor- 
rection equation, based on the SIMPLE-C algorithm. 14 
In constructing cell face fluxes, a momentum interpola- 
tion scheme 15 is employed which introduces damping in 
the continuity equation. Grid motion/deformation terms 
are implemented in a Geometric- Conservation-Law (GCL) 
preserving fashion. 16 A dual-time formulation is employed 
where at each physical time step, between 5 and 20 pseudo- 
timesteps of the SIMPLE-C algorithm are applied. Specif- 
ically, at each pseudo-timestep, the discrete momentum 
equations are solved approximately (using a simple point 
iterative scheme), followed by a more exact solution of the 
pressure correction equation (using the PETSC 17 parallel 
LU preconditioning+GMRES utilities). Turbulence scalar 
and energy equations are then solved in succession. Paral- 
lelization is implemented in a standard fashion by invoking 
domain decomposition based on METIS 18 in the front end, 
and MPI based message passing in the CFD code. All 
of the large scale simulations presented in this paper were 
executed on the Columbia supercomputer at NASA Ames 
Research Center. The code scales very well on this system 
as illustrated in Figure 2. 



Fig. 2 Parallel efficiency of NPHASE-PSU on Columbia sys- 
tem for a l.lxl 0 6 cell spinning cylinder case. 


DiRTlib and SUGGAR 

The overset grid approach 19 utilizes a composite grid 
consisting of a set of overlapping component grids to 
discretize the domain. No point-to-point or face-to-face 
matching is required between component grids. The so- 
lution on the component grids are linked by identifying 
appropriate intergrid boundary points (IGBPs) where the 
solution is given by a specified boundary value obtained 
by interpolation from another overlapping donor compo- 
nent grid. Tire overset domain connectivity information 
(DCI), which consists of the identification of the intergrid 
boundary points and corresponding interpolation sources, 
is obtained by an overset grid assembly step. The current 
effort utilizes two overset software libraries to add the over- 
set capability to NPHASE-PSU. 

DiRTlib, 20 which stands for Donor interpolation Recep- 
tor Transaction library, is a solver neutral library that en- 
capsulates the functionality required by the solver to utilize 
the overset domain connectivity information. It is inde- 
pendent of the solver grid storage and topology, dependent 
variables, etc. and can be used with any solver. The solver 
calls a few functions to initialize the library, load the DCI 
interpolation, and transfer the data to appropriate proces- 
sors in a parallel execution environment, and apply the in- 
terpolated data as boundary values at IGBPs. Solver func- 
tions must be provided and are called by DiRTlib to get and 
put data in the correct solver dependent variable storage 
locations. When the solver executes in a distributed mem- 
ory parallel computational environment the solver must 
also inform DiRTlib of the parallel decomposition enabling 
DiRTlib to get/put data from the appropriate parallel pro- 
cess. 

The current overset grid assembly process is performed 
using the SUGGAR code, 21 which stands for Structured, 
Unstructured, Generalized overset Grid AssembleR. SUG- 
GAR is a general overset grid assembly code with the ca- 
pability to create the domain connectivity information at 
node and/or element centers for most current grid topolo- 
gies including any combination of structured Cartesian and 
curvilinear, unstructured tetrahedral and mixed element, 
general polyhedral, and octree based Cartesian grids. For 
static grid assemblies with no motion between component 
grids the grid assembly is a pre-processing step. The case 
of solution and time dependent motion, requires the solver 
to communicate the new body and grid positions to the grid 
assembly process, wait for it to complete, and then load the 
new DCI. For the case of prescribed motion, such as used 
in the present study, the DCI is a-priori computed and saved 
in a file for each time step in the simulation and the solver 
simply loads the file appropriate for each time step. 

The donor interpolations produced by SUGGAR are 
a set of linear weights that multiply the values at the 
donor members. For a cell centered flow solver, such as 
NPHASE-PSU, the interpolation stencil will use as mem- 
bers the cell in the donor grid that was found to contain an 
IGBP and the neighboring cells that share a face with tire 
donor cell. The interpolation weights are computed using 
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an unweighted least square procedure. 


Results 


Verification Studies 

In the context of the overset meshing strategy employed 
for gearbox windage simulations, the meshes will be in 
motion relative to one another. As indicated above, the ap- 
proach taken here is to solve the flow in the absolute frame 
of reference for the entire computational domain, i.e. on 
all meshes, those that are rotating and those that are sta- 
tionary. In order to verify that NPHASE-PSU correctly 
handles these gear relevant rotating mesh systems, two ver- 
ification studies were performed, rotating Couette flow, and 
flow near an infinite rotating disc, both of which have avail- 
able analytical solutions. 

Rotating Coutte Flow 

Figure 3 is an illustration of the 33x81 (radial x az- 
imuthal) computational domain for the incompressible ro- 
tating Couette flow case. The inner and outer boundaries 
are walls. In this case, rimer = 0-5 and r 0U ter= 10. The 
outer cylinder is held stationary and a rotation rate of 
© = Is -1 is specified for the inner cylinder, The Reynolds 
number independent analytical solutions for the tangential 
velocity and pressure are: 




P 

C 




(0 

3 


(4) 


(5) 

( 6 ) 


Figure 4 shows a comparison of the analytical solution with 
three NPHASE-PSU runs, designated Approaches 1, 2, and 
3. Approach 1 solves the absolute velocities in the ab- 
solute frame on a stationary grid (adapting inner cylinder 
boundary conditions accordingly). Approach 2 solves for 
the relative velocities in the relative frame on a stationary 
grid (i.e. frame of reference rotating with angular- velocity, 
or, adapting the momentum equation source terms and outer 
cylinder boundary conditions accordingly). Approach 3 
solves for absolute velocities in the absolute frame using 
a time accurate analysis on a rotating mesh. Approach 3 is 
the most relevant for gear analysis. Figure 4 illustrates that 
the code returns the analytical solution for all three simula- 
tion approaches to within the accuracy of the second order 
accurate discretization numerics and grid used. 

Flow Near an Infinite Rotating Disc 

The second verification case is a classic 3-D exact solu- 
tion to the incompressible Navier- Stokes (N-S) equations. 
An infinite radial span disk rotates with an angular velocity, 
©. This induces tangential flow in the direction of rotation, 
radial outflow, and an axial flow towards the center of the 
disc. In this case, the N-S equations reduce to a system of 



Fig. 3 Computational domain for the rotating Coutte flow 
case. 



Fig. 4 Comparison of NPHASE-PSU and analytical solutions 
for the rotating Couette flow case. 


non-linear ODEs which can essentially be exactly solved 
numerically. Figure 5 is notional sketch of the flow field 
from Schlichting. 22 

Figure 6 shows the 232662 element unstructured mesh 
employed for the analysis. The extent of the computational 
domain was selected to be r max = 1.0, Zmax = 1.0. For a 
choice of co = 1.0 and j u = l.OxlO -2 , this domain provided 
that the solution sampled within the region r < 0.2, z < 0.4 
compares very closely with the exact solution despite the 
necessarily finite extent of the domain, 

NPHASE-PSU was applied using two approaches. Ap- 
proach 1 solves the absolute velocities in the absolute frame 
on a stationary grid (adapting disk boundary' conditions ac- 
cordingly). Approach 2 solves for absolute velocities in 
the absolute frame using a time accurate analysis on a ro- 
tating mesh. Figure 7 illustrates that the code returns the 
exact solution for both simulations to within the accuracy 
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Fig. 5 Sketch of the flow field in the vicinity of an infinite span 
rotating disc. 



Fig. 6 232662 element unstructured mesh employed for the 
infinite spanning rotating disc. 


of the second order accurate discretization numerics and 
grid used. 

Couette Flow with Wall Heating 

The third validation performed sought to verify the 
viscous dissipation term implementation in NPHASE- 
PSU. The relevance of viscous dissipation physics to gear 
windage is significant as discussed above. The analytical 
solution for the temperature distribution in a laminar linear 
Couette flow was chosen. Figure 8 shows a diagram of the 
configuration. The product of the Prandtl number (TV = 
fiCp/k) and the Eckert number (Ec = Uf /C P (T\ - 7b)), 
PrEc , is a measure of the role of viscous dissipation in 
a flow. The nearly exact comparisons between CFD re- 
sult and the analytical solution, shown in Figure 9, across a 
range of PrEc, illustrates that the viscous dissipation terms 
in NPHASE-PSU are implemented correctly. 

Validation Studies 

The experimental data of Diab et al. 8 was used to val- 
idate NPHASE-PSU for the case of unshrouded, isolated, 
rotating spur gears. Diab et al. s tested four different spur 
gears and a disk in free air on a spin-down test rig. The 
gears varied in diameter, width, and tooth count. Tire prop- 


• F = V /ro. analytical solution 
G e Mjra. analytical solution 

• -H - -V I /vy1|\t,)), analytical solution 
F - V> M, NPHASE-PSU, stationary mesh 

G = MJm, NPHASE -PSU, stationary mesh 


1 - -H - -VJsrj gwa). NPHASE-PSU, stationary mesh 

F = V,/ru, NPHASE-PSU, rotating meah 

- \ G a \IJra. NPHASE-PSU. rotathg mesh 



Fig. 7 Comparison of NPFIASE-PSU and exact solutions for 
the infinite span rotating disc. 


U=U„ T=T| 



11=0, T=T 0 

Fig. 8 Illustration of the linear heated Couette flow case. 



Fig. 9 Comparison of NPHASE-PSU and exact solutions for 
the linear heated Couette flow case. 


erties of the gears and disk used by Diab et al. 8 are provided 
in Table 1 . Diab et al. 8 did not study the effects of gear en- 
closure or lubrication. A sequence of prescribed constant 
rotation rate simulations was used to replicate the experi- 
ment. 
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Table 1 Diabet al. 8 Gear Properties. 



Pitch Diameter 
(mm) 

Tooth Width 
(mm) 

Module 

(mm) 

Gear 1 

288 

30 

4 

Gear 2 

144 

30 

4 

Gear 3 

144 

60 

4 

Gear 4 

144 

60 

6 

Disk 

300 

30 



number rotation speeds. All cases were run for at least 2 
complete revolutions to remove simulation startup transient 
behavior. Convergence histories show that transients leave 
the solution after about one revolution as illustrated in Fig- 
ure 11, where it is also observed that pseudo-time residual 
drop approximately 2 orders of magnitude in each physical 
timestep when 10 pseudo-timesteps are used per physical 
timestep. 


10- F 


Single Grid Simulations 

Grids were generated for all four spur gears and the disk. 
For the gear studies where the high Reynolds number k- 
e turbulence model was used, near- wall grid spacing was 
defined to accommodate wall-functions (e.g., y + ~ 70 for 
gear 1, <o= 1000 .? _1 ). The single plane of symmetry in the 
problem was exploited to reduce total cell count by a factor 
of 2. Grid cell counts for the different cases varied between 
2x 10 6 (Gear 4) and 8xl0 6 (Gear 1). Grid generation was 
further simplified by employing a hybrid mesh topology, as 
illustrated in Figure 10. Specifically, for the regions above 
the surface of the gear teeth to the outer boundary, struc- 
tured hexahedral cells were used. For the region above the 
gear face surface, unstructured prism cells were employed. 
The meshes were generated using the commercial grid gen- 
eration software package Gridgen. 23 The computational 
domain of the isolated gear grids was extended to approx- 
imately five times the gear radius from the gear surface in 
all directions. It was found that this distance was adequate 
for defining a symmetry boundary condition since the flow 
is nearly stagnant there. 

An azimuthal step size of 1 /40th of one tooth passage 
duration (the time it takes one tooth to rotate to the position 
of the tooth adjacent to it) was used in all CFD calculations. 
This corresponds to 2880 timesteps per gear revolution for 
Gear 1, 1440 timesteps for Gears 2 & 3, and 960 timesteps 
for Gear 4. All cases used 10 pseudo-time iterations per 
physical timestep. 

CFD runs were made for four gears and the disk at a 



Fig. 10 Grid topology of Gear 4. 



iteration 


10 pseudo-timesteps 
: per physical timestep 



irv — 1 i i i 1 i i I i i i i I i— -i i i I 

lyOOO 19050 19100 19150 19200 

iteration 


Fig. 11 Example convergence history for Diab Gear 4. 

Comparisons between the power loss results of Diab 
et al. 8 and the NPHASE-PSU analysis are presented in 
Fig. 12 for all four gears. The CFD analysis for all 4 gears 
exhibited very good agreement with experiment. The disk 
case, however, did not share this same level of agreement, 
as illustrated in Figure 13, where NPHASE-PSU results 
are seen to underpredict the measured power loss. In or- 
der to elucidate the reasons for the deterioration in solution 
accuracy observed for the disk case, a number of observa- 
tions and studies were made. First, it is observed that the 
measured (and computed) windage loss power levels for 
the disk are much smaller than the comparably sized spur 
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gear (Gear 1, D ~ 300 mm). This arises due to the ab- 
sence of any azimuthal pressure variation in the disc flow - 
torque losses are due entirely to viscous effects, and these 
are clearly underpredicted. Indeed the absolute magnitudes 
of loss underprediction between the disk and Gear 1 are 
comparable (e.g. ~ 50W @ 600s _1 ), so presumably this 
underprediction of shear is present in all of the gear simu- 
lations, however its relative magnitude is small. 

To explore this further, the low-Reynolds number Menter 
model was applied to the spinning disc case (using an ap- 
propriate sublayer resolved mesh). Figure 13 illustrates 
that improved turbulence modeling does benefit solution 
accuracy especially at higher rotation rates. This observa- 
tion is not as important for the gear cases since the spin- 
down torques associated with azimutally varying pressure 
forces in the vicinity of the gear teeth dominate the viscous 
forces as shown in Figures 14 and 15. It can be seen that, 
especially at lower rotation rates, the contribution of vis- 
cous loss to total windage loss, is small. However, we do 
observe that at higher pitchline velocities tire relative mag- 
nitude of the viscous torque increases with respect to the 
pressure torque. This can be see for higher rotation rates 
on Gear 1 (Figure 14) as well as by considering the smaller 
size (and hence lower pitchline velocities) of Gear 2 (Fig- 
ure 15). 

In summary our results suggest that for the very high 
speed gears to be encountered in rotorcraft (and other high 
performance aircraft) transmissions, viscous effects will 
become more important than encountered here, and will 
therefore require attendant research attention. 



Fig. 12 Comparisons between the experimental results and 
the NPHASE-PSU analysis. 

Details of 3-D Flow Field 

A number of important physical features of the predicted 
flow field are available upon interrogation of the CFD simu- 
lations. Figure 16 shows a view of the predicted secondary 
velocity vectors on the symmetry plane in the gear relative 
frame of reference for one of the Diab cases. One can see 



Fig. 13 Effect of turbulence mode) selection on viscous work 
prediction. 
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Fig. 14 Breakdown of windage pow er losses for Diab 1. 


a significant vortical structure within the gear tooth region, 
and the tooth-to-tooth periodicity that has been achieved 
in the transient simulation. Figure 17 shows a view of the 
predicted surface pressure distributions for the four Diab 
gears. There, comparatively large vs. small pressures are 
observed on the leading and trailing tooth faces - this dif- 
ference being the source of tire pressure component of the 
spin-down torque. The highly three-dimensional nature of 
the flow in these relatively low aspect ratio spurs gears are 
also clearly seen (figure shows only 1/2 of each gear). Sig- 
nificant 3-D effects are also clearly visualized in Figure 18 
where gear relative streamlines are displayed in the near 
tooth region along with an isosurface of predicted static 
pressure in a region of high pressure. 
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Fig. 15 Breakdown of windage power losses for Diab 2. 



Fig. 16 Predicted velocity vectors in the symmetry plane of a 
Diab gear. 



Fig. 17 Predicted surface pressure distributions for the four 
Diab gears. 


Overset Grid Simulations 

As of the writing of this paper, the overset capabil- 
ity for spuming gear simulations has been established in 
NPHASB-PSU. In this context, we are pursuing an overset 
verification effort and a relevant validation effort. The veri- 
fication effort involves solving isolated spinning gear cases 
studied above using a rotating near-gear mesh and a station- 
ary farheld mesh, the necessary approach to be used in on- 
going gearbox windage activities. Figure 19 shows a view 
of an overset Diab Gear 4 simulation. Surface pressure 



Fig. 18 Predicted relative frame streamlines and isosurface 
of a high pressure region for Diab Gear 4. 


contours are displayed along with overset gear and back- 
ground meshes on the symmetry plane. As expected, the 
CFD code returns nearly identical results for this moving 
grid overset simulation to the non-overset results reported 
above. We continue to parameterize gridding requirements 
in the overlap region, including the proximity of the over- 
lap region to the gear in assessing the retained accuracy of 
the overset approach. 



Fig. 19 Verification study: 3-D overset grid solution and 
topology for Diab Gear 4. 


The validation effort under way involves simulating a se- 
ries of shrouded low speed gears for which windage loss 
measurements were made by Dawson. 7 To date we have 
developed valid overset assembly grid topologies for this, 
as shown in Figure 20, and anticipate reporting CFD anal- 
ysis for these configurations in the near future. 


Conclusion 

This paper has summarized the adaptation and validation 
status of a CFD method for gear windage aerodynamics. 
Validation studies of 3-D spur gears in free space demon- 
strate very good agreement with published data. The fol- 
lowing conclusions apply: 1) Viscous/turbulence modeling 
was identified as a shortcoming since loss power was con- 
sistently underpredicted for the viscous-drag-only spinning 
disk cases. 2) Low Reynolds number (sublayer resolved) 
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Fig. 20 Validation study: 3-D overset grid topology for Daw- 
son shrouded gear. 


turbulence modeling exhibited improved performance for 
this case. 3) High Reynolds number two equation mod- 
eling proved appropriate for modeling the moderate speed 
Diab spur gear suite, due apparently to the dominance of 
the pressure torques on spin -down. 4) The budgets of vis- 
cous and pressure components of spin-down torque suggest 
that viscous effects will bee ome much more important, per- 
haps exceeding 50 percent, for gears with pitchline speeds 
approaching twice that investigated here. 5) Overset mesh- 
ing is and will continue to be a critical enabler in this effort. 
The capability has been established and demonstrated here 
and will become integral in further studies where the iso- 
lated gear assumption will not be relevant. 

Our continuing work in the gear windage aerodynam- 
ics area will focus on: 1) Further validation for shrouded 
gears, 2) Turbulence modeling including sublayer mod- 
eling, v 2 -f modeling and/or Reynolds Stress Modeling 
(RSM) and/or Detached Eddy Simulation, 3) The role of 
laminar-turbulent transition, 4) Multi-phase flows, 5) Over- 
set gridding in the context of gear tooth contact, 6) Indus- 
trial and NASA gearbox configuration applications, and 
7) Evolution of design guidance for minimizing windage 
losses 
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An important and challenging element of CFD modeling for gear systems is the relative motion 
between, and the contact between individual gears. This challenge led us to invent/develop what we 
have now come to term the Hybrid Overset Immersed Boundary Method (HOIBM). The HOIBM was 
developed under this NRA program and was first published by us at the AIAA CFD Conference in 2009. 
Although that paper, AIAA-2009-3992, included many other recent overset grid technology 
advancements, we include below only the content of the paper that treated the HOIBM. In addition to 
the technical explanation and canonical test case presented in that paper, we also presented to the NASA 
sponsor, and VLC and VLRCOE audiences, the application of the method to spinning spur gears in 
contact. Some of these results are included here (although they can be best appreciated only by viewing 
an animation of the time accurate simulation of the gear meshing). 



Elements of HOIBM application to meshing spur gears. Grids, overset assembly and solution time 

snapshots. 
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(Extracted from AIAA-2009-3992): 


VIII. Demonstration of Immersed Boundary Approach to Eliminating 

Orphans 

The use of the Suggar+-t- identification of immersed boundary locations to eliminate orphans in an overset 
grid system will be demonstrated for two different problems using two different flow solvers. The results 
presented here are intended as demonstrations only and additional development and validation of the results 
will be presented in the future. 

A. Two-Dimensional Glottal 

The first demonstration problem is a simple two-dimensional geometry. This configuration, which is illus- 
trated in Figure 47, is a canonical model of the human glottis, for which experimental measurements are 
available. 19 In the experiment, the nominal “vocal folds" are closed at t=0. and a pressure difference is 
applied between the inflow and outflow boundaries. The “vocal folds” are then opened and closed over a 
period of 5.G7 seconds, following a non-sinusoidal path controlled by a stepper motor. As the folds open, the 
Ap across the test section drives a jet flow through the space between the folds. This test case challenges the 
traditional overset approach due to the point contact between the cylindrical folds, which would normally 
result in orphans. This motivated the mesh refinement near the contact point illustrated in the figure, so as) 
to suitably approximate point closure. 

The equations for the fluid motion was solved using the NPHASE-PSU code, an unstructured, parallel 
finite volume solver. Algorithmically, NPHASE-PSU follows the well-established segregated pressure based 
method. A colocated variable arrangement is used, and a lagged coefficient linearization is applied. In this 
laminar flow example, 2 nd order upwind and central discretization schemes were selected for the convection 
and viscous terms in the momentum equations. Continuity is introduced through a pressure correction 
equation based on the SIMPLE-C algorithm. In constructing element face fluxes, a Rhie-Chow momentum 
interpolation scheme is employed. Further details of the scheme are available in Reference 20. DirtLib and 
the necessary “iblanking” instrumentation for the hybrid overset- immersed method w'ere recently installed in 
the code. Specifically, all out-immersed elements are nullified through appropriate source term treatments, 
and wall faces are inserted in the data structure at all internal faces that separate out-immersed from active 
elements. 

Since the glottal motion was specified, all of the overset assembly files were pre-generated using SUGGAR 
driven by a Python motion controller. In order to accommodate the fully closed conditions (no flow, despite 
a driving Ap between inlet and outlet), a constant stagnation pressure condition was specified at the inlet to 
match the experimentally-obtained peak jet centerline velocity. A constant reference pressure was specified 
at, the domain outlet for this incompressible flow. A dual-time approach was used w'herein 50 sub-iterations 
were employed for each physical time step. At each sub-iteration, the discrete momentum equations are 
solved approximately, followed by a more exact solution of the pressure correction equation. 

Figure 47 shows a portion of the overset grid system for this case when the Glottal is partially open. 
The grid for the upper glottal is displayed in red, the channel grid is shown in gray, and the tw T o lower collar 
grids are showm in cyan. The boundaries of all overset grids are also displayed with thick fines. Figure 48 
shows the complete overset grid system w T hen the gap is closed and the surfaces are in contact. 

Figure 49 shows a closeup of the grid system near the contact point. The upper and lower Glottal grids 
are shown in red and blue respectively. The cyan and magenta dots are located at the element centers of 
fringe elements. The red and blue dots are located at the center of elements marked as OUT-IMMERSED 
in the upper and lower Glottal grids respectively. The faces between active and OUT-IMMERSED elements 
are showm as thick lines and indicate the approximate geometry being used in the solution at this time step. 

Figure 50 shows a sequence of the predicted instantaneous axial velocity through one open-close cycle. 
The experimental configuration has a square cross section at the inlet, rather than an infinite spanwise 
extent assumed with the present 2D simulation. The authors are presently extending the model to capture 
the three-dimensionality of the geometry and the flow field. 
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Figure 47. Portion of overset grid system for the 2D Glottal case. 



Figure 48. Complete overnet grid system for the 2D Glottal case when the gap in cloned. 



Figure 49. Closeup grid system near the contact point for the 2D Glottal case when the gap is closed. 
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As the first phase of the NRA program wound down and we were funded by NASA to pursue a 
Phase 2, we turned our attention to further validation, physics interpretation/understanding and design 
guidance. We had for the first time available to us the recently obtained NASA Glenn Windage Facility 
data for validative comparison. We moved to another CFD code, OVER-REL, for this work due to its 
well established efficiency for single tooth periodic boundary condition capability. Parametric shroud 
configuration studies carried out in the Glenn experiments and the CFD analyses elucidated important 
physical mechanisms of windage losses as well as suggested mitigation strategies due to shrouding and 
newly proposed tooth contour modifications. The NASA Glenn data and CFD analyses showed that 
axial and radial gear shrouding are effective in significantly reducing gear windage losses both 
independently and when employed together. For all configurations studied, the dominant physical 
mechanism contributing to windage losses is the pressure field associated with diversion and 
impingement of the high speed relative flow on the leading tooth surface. This was a discovery not 
understood by the community prior to our research, ft was found that shrouding mitigates the magnitude 
of this mechanism but not its dominance in the loss budget. Although viscous effects were identified as a 
secondary source of windage losses, they are not insignificant. Also, the sublayer resolved turbulence 
model used, although more accurate than a high Reynolds number form with wall-functions, apparently 
still underpredicts viscous losses for the rotational boundary layers studied. Viscous losses were found 
to be reduced by shrouding as well. However, small axial clearances were observed computationally to 
increase viscous losses compared with larger shroud clearances, suggesting that this effect could become 
important for even smaller clearances and at higher speeds. The CFD results showed good agreement 
with the Glenn experiments which exhibited quite similar loss trends to the idealized shrouded gear 
configurations studied computationally. T 

The physics understanding that our CFD studies led to, suggested a set of possible geometric 
tooth modifications to further reduce windage loss. We studied these notional concepts and one of them, 
a spoiler or ramp protruberance appeared quite promising. Further study of this concept led to our 
patenting the concept. The 2010 AHS Forum paper summarized this early research. We submitted the 
paper to the ASME Journal of Fluids Engineering and it was accepted and appeared in March 2011. This 
JFE version is included here in its entirety. 
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CFD Analysis of Gear Windage 
Losses: Validation and 
Parametric Aerodynamic Studies 

A computational fluid dynamics (CFD) method has been applied to gear configurations 
with and without shrouding. The goals of this work have been to validate the numerical 
and modeling approaches used for these applications and to develop physical under- 
standing of the aerodynamics of gear windage loss. Several spur gear geometries are 
considered, for which experimental data are available. Various canonical shrouding con- 
figurations and free spinning (no shroud) cases are studied. Comparisons are made with 
experimental data from open literature, and data recently obtained in the NASA Glenn 
Research Center Gear Windage Test Facility, Cleveland, OH. The results show good 
agreement with the experiment. The parametric shroud configuration studies carried out 
in the Glenn experiments and the CFD analyses elucidate the physical mechanisms of 
windage losses as well as mitigation strategies due to shrouding and newly proposed 
tooth contour modifications. [DOI: 10.1115/1.4003681] 


1 Introduction 

Gear windage losses refer to the power losses that are a result 
of aerodynamic forces (pressure and viscous) acting against the 
rotation of high speed gears. Windage losses are a source of sig- 
nificant heating and fuel consumption in rotorcraft and other ver- 
tical take-off and landing (VTOL) systems. The weight and pack- 
aging constraints inherent in these systems require the gearing 
components to be both lightweight and heavily loaded. Attendant 
to this, the gears are required to operate at high rotational speeds, 
where windage losses become significant with respect to other 
gearbox losses (rolling, sliding, and lubrication losses). 

A host of experimental studies of gear windage have appeared 
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in literature [1-10], These studies employ either isolated gears 
[1-5] or closed-loop gear-train systems [6-10], Generally, these 
studies parametrize gear geometry elements, rotational speed, en- 
closure geometry, and lubrication system characteristics (flow 
rate, jet location, and lubricant rheology), and use dimensional 
analysis to develop correlations for the power losses that have 
proven useful in the design process. 

It was shown experimentally by Dawson [1 ,2] and Winfree [3] 
that geometric modifications of the near-gear flow path, such as 
with shrouding, can significantly reduce both windage losses and 
lubricating oil consumption (80% and 40% reductions observed, 
respectively, in Ref. [3]). Such modifications are of significant 
interest in rotorcraft transmission design. Unfortunately, only a 
few empirical studies have been made available to date. Accord- 
ingly, there is a need for more systematic and prototypical experi- 
mental data and a need for improved understanding of the physi- 
cal mechanisms involved in the schemes used to reduce these 
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aerodynamic losses. 

Recently, the NASA Glenn Research Center has installed a gear 
windage test facility to address these issues [5], This facility is 
expected to produce the most comprehensive set of gearbox wind- 
age data to date, and is providing CFD validation data required in 
the present NASA sponsored research effort. The test facility is 
designed to parametrize the effects of spur gear geometry, shroud 
geometries, different lubrication system configurations, system 
pressures and temperatures, and gear meshing on windage loss. 
Data from both open literature data and new NASA Glenn facility 
are used in this paper. 

CFD practitioners have also recently begun to analyze gear 
windage aerodynamics toward the goal of developing a more gen- 
eral design approach. Recent 2D studies [11] were performed us- 
ing the commercial flow solver fluent, where a side correlation 
factor was used to account for 3D effects (although these authors 
stated that work is under way to extend their simulations to 3D 
[12]). Marchesse et al. [13] used a similar approach as Al-Shibl et 
al. [11] but used a structured grid with the flow solver ansys cfx. 
They also extended their 2D model into three dimensions, using 
one tooth passage with periodic boundary conditions as their com- 
putational domain. Imai et al. [14] investigated 3D bevel gears in 
mesh in an air-oil atmosphere, modeling the gears as porous bod- 
ies. Hill et al. [15] implemented a 3D, unstructured, overset mov- 
ing mesh method and applied it to isolated spur gears in air and 
validated the method against data of Diab et al. [4], Most recently, 
Hill [16] applied a 3D, structured, overset, steady state CFD 
method to elucidate the physical mechanics of how shrouding 
reduces windage effects on spur gears. 

The validated capability developed in our earlier work [15] has 
motivated, among other things, the focus of the present study: to 
use a CFD model to interrogate the physics associated with gear 
windage aerodynamics. The attendant goal of this is to develop 
design guidance principles for methods to minimize these power 
losses in rotorcraft, and, potentially, other high speed gear sys- 
tems. This has been and will continue to be enabled by the avail- 
ability of the NASA Glenn Research Center data sets. 

This paper is organized as follows. First, the governing equa- 
tions and numerical procedures are reported. A series of numerical 
studies are then performed, which provide insight into several of 
the physical mechanisms of windage losses, including the role of 
gear shrouding. This is followed by an overview of the NASA 
Glenn Gear Windage Test Facility with a selection of data ob- 
tained to date. CFD simulations of the Glenn windage loss tests 
are subsequently presented and discussion of these comparisons is 
provided. Finally, the results of numerical experiments, guided by 
the findings described in this paper, are presented, which demon- 
strate a benefit in windage performance beyond that observed with 
shrouds alone. 


2 Theoretical Formulation 

2.1 Governing Equations. This paper is focused on the 
physical mechanisms of gear windage loss, and mitigation 
schemes for reducing them. The available data [1-5] show that 
windage losses already become significant at low subsonic tip 
Mach numbers. Accordingly, in this work, we restrict ourselves to 
incompressible analysis. Also, as will be shown, much can be 
learned from systems that are either shroud-free or have com- 
pletely axisymmetric shrouds (i.e., fully enclosed radial shrouds), 
so we consider single-tooth domains that are periodic and steady 
in the frame of reference of the spinning gear. The conservation of 
mass and momentum can be written in integral conservation law 
form for a system rotating with constant angular velocity ® as 

pV-dS = 0 (1) 

J s 


P VV-dS = - pdS + T-dS - pa>X (coX r)d V 

Js Js Js Jv 


■1 


2 p{w X V)dV 


(2) 


where V is the velocity vector in the relative frame of reference. It 
was shown in Ref. [15] that a sublayer resolved two-equation 
turbulence model performed better than a high Reynolds number 
form with wall-functions in predicting viscous losses on a spin- 
ning disk. Accordingly, in this work, we adopt T=(p,+ 

+ (Vy) r ) and a sublayer resolved q-co turbulence model due to 
Coakley [17] is used. The dependent variables in this model are 
related to the turbulence kinetic energy k and the turbulence dis- 
sipation rate s through q='Jk and co=s/k. In this model, the eddy 
viscosity is obtained from 

ft = pC t ,Dq 2 /a) (3) 

where C^= 0.09 and D is the near-wall damping function 

D = 1 - e-owVt* (4) 

where a= 0.02 and d„ is a measure of the normal distance to the 
nearest wall. The modeled transport equations for q and a> are 

(5) 

J po)Y ■ dS = J + “j V (O ■ dS + J piQC^S - C 2 <o 2 )dV 

(6) 


2.2 CFD Code and Numerics 

2.2.1 over-rel. over-rel is a parallel, cell-centered, incom- 
pressible, finite-volume code based on structured overset multi- 
block grids and the time-marching, pseudocompressibility formu- 
lation of Chorin [18], Inviscid fluxes are formulated from the 
Roe-approximate Riemann solver [19] and extended to third-order 
accuracy through the monotone upstream-centered scheme for 
conservation laws (MUSCL) scheme [20], Second-order accurate 
central differences are utilized for the viscous fluxes. Numerical 
derivatives are used to calculate the flux Jacobians. A symmetric 
Gauss-Seidel method is applied to solve the resulting linear sys- 
tem of equations. In the present application, the code’s turboma- 
chinery analysis instrumentation is employed; all simulations are 
carried out for a single-gear tooth, with periodic boundary condi- 
tions, in a noninertial frame of reference rotating with the gear. 

2.2.2 Overset Meshing: sugoar and dirtlib. The tight shroud 
clearances studied lead to mesh topology constraints that, in turn, 
lead to poor quality block-structured meshes, over-rel has over- 
set meshing capability, which enables high quality meshes for 
these configurations. The overset assembly process is performed 
using the suggar code [21], a general overset grid assembly code 
with the capability to create the domain connectivity information 
(DCI) at node and/or element centers for general grid topologies. 
For static grid assemblies with no motion between component 
grids (as here), the grid assembly is a preprocessing step. The 
intergrid interpolations produced by suggar use an unweighted 
least square procedure. 

over-rel is instrumented with dirtub [22], a solver neutral 
library that encapsulates the functionality required by the solver to 
utilize the overset DCI produced by suggar. dirtlib is indepen- 
dent of the CFD solver’s data structures, and can therefore be used 
with any solver. The CFD code calls a few functions to initialize 
the library, load the DCI interpolation, transfer the data to appro- 
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Table 1 Ref. [4], gear properties 



Pitch diameter 
(mm) 

Tooth width 
1 nun I 

Module 

(mm) 

Teeth 

Gear 1 

28S 

30 

4 

7 1 

Gear 2 

144 

30 

4 

36 

Gear 3 

144 

60 

4 

36 

Gear 4 

144 

60 

6 

24 

Disk 

300 

30 




prude processors in a parallel execution environmem. and apply 
the interpolated data as boundary values at intersrid boundary 
points. 

3 Results 

3.1 Loss Mechanisms and lludgets in Gear Windage. In 
their 21X14 papet. [>iab el al. [4] studied four spur gears and a disk 
spinning lively in air (i.e.. no shrouding). The gears varied in 
diameter, w idth, and tooth module. Table l summarizes the gear 
properties. In the present work, free spinning simulations were 
first carried out for Diab gear I and the disk using ovitK-Kitt.. 
Nono verse! multi block-slntelii red meshes were employed for 
both, as shown in Fig, I Near-Wall grids were constructed to 
return wall cell y* values <1 and wall normal stretching ratios 
<1.2 everywhere in order to adequately resolve the high Rey- 
nolds number boundary layers that arise. The grid topologies, 
near-wall and spacing, and grid stretching ratios were maintained 
as closely as possible between the gear and disk meshes. Per the 
periodic boundary conditions employed, one tooth passage 
(2rr/72) is modeled for both configurations In order to stably 
time-march the overrsl solution, a very small inflow velocity 
and a pressure outflow boundary were included adjacent to the 
maximum and minimum axial boundaries upstream and down- 
stream of the miming elements. This artificial through-flow veloc- 
ity was successively reduced to where no perceptible change in 
loss values was relumed. 

Figure 2 shows ovrat-Riit. results for gear I with experimental 
data and previous CPU results by the authors [15]. The OVF.R-RE1. 
results show very good agreement with both experiment aud pre- 
vious simulations Figure 3 shows that for gear 1. the pressure 
torque associated with the integrated pressure difference between 
leading aud trailing loolli surfaces dominates the loss budget As 
the rotation rale increases, viscous losses remain a nearly constant 
fraction of total loss (1091'). Figure 4 show s a comparison of loss 
results between gear 1 and the disk, lhe total losses arc much 
smaller lor the disk and ;ue due to viscous shear alone. 

In Fig. 5. lhe torque per unit span contributed by viscous shear 
is compared for the disk (up to its outer radius) and gear 1 (up to 
its base radius), file geometry of these systems requires that the 
all pressure spin-down torques are due lo the pressure differences 
between the leading and Railing tooth surfaces. Figure 5 illustrates 



Fig. 1 Comparison of Diab gear 1 and disk grids 

Journal of Fluids Engineering 



Fig. 2 Comparison of results from experiment (.4.1, NPHASE- 
PSU [15], and over-rel for Diab gears 1-4 


dial the viscous losses are very similar between these configura- 
tions (indeed, tile small differences in predicted viscous power 
values are due to very small grid differences) showing, conclu 
sively, that 3D effects (i.e., nonaxisymmetric) associated with 



Fig. 3 Comparison of results from experiment [4] and over-rel 
for Diab gear 1. Viscous and pressure loss budgets are 
included. 



Fig. 4 Comparison of Diab gear 1 and disk measurements and 
over-rel solutions 
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Fig. 5 Comparison of predicted torque per unit span contrib- 
uted by viscous shear for the Diab disk (up to its outer radius) 
and gear 1 (up to its base radius) 


pressure forces are almost completely responsible for the signifi- 
cantly greater losses for the gear. 

The physical mechanisms associated with the dominant pres- 
sure torque can be studied by interrogating the CFD results for 
gear 1 . In Figs. 6-8. several 3D visualizations are presented for 
the 850 rad/s simulation. In Fig. 6. a number of relative frame of 
reference streamlines, colored by local static pressure, are dis- 
played. These streamlines are seeded close to the gear face and 
teeth and integrated in both directions. Some of the high speed (in 
the relative frame) tangential flow near the gear face plane is 
diverted into the tooth passage, where strong secondary flows are 
evident. A stagnation region is observed on the leading tooth sur- 
face. where the diverted relative flow impinges near the gear face. 



Fig. 6 Three-dimensional relative frame streamlines colored 
by static pressure for Diab gear 1 at 850 rad/s 



Fig. 7 Axial projection of velocity vectors halfway between 
gear face and gear centerline for Diab gear 1 at 850 rad/s. Vec- 
tor density of 0.5. Background contours of local normalized 
projected relative velocity magnitude. 
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Fig. 8 Axial projection of velocity vectors near-gear centerline 
for Diab gear 1 at 850 rad/s. Vector density of 0.5. Background 
contours of local normalized projected relative velocity 
magnitude. 


By symmetry, this axial transport arises on both sides of the gear 
and therefore leads to impingement of oppositely directed flow 
and radial ejection of momentum near the gear centerline. An 
axial projection of relative velocity vectors is shown in Fig. 7 at a 
plane halfway between the gear' face and the gear centerline. A 
vector density of 0.5 (vector plotted at every other grid point) is 
applied for clarity. Contours of local no rmalize d projected relative 
velocity magnitude are included (V' : =\IV*+V~/ wr). Two counter- 
rotating passage vortices are present. Peak normalized secondary 
velocity magnitudes near l/2«r are observed, indicating the 
strength of these secondary motions. The flow at this axial loca- 
tion is reminiscent of a rearward facing step and/or cavity flow 
with attendant vortical recirculation regions. Figure 8 shows the 
same plot but near the gear centerline. Here, one sees the very 
significant radial ejection quite clearly. The flow has a component 
directed upstream (against the relative flow) near the leading sur- 
face at the tip radius. The magnitude of the ejection flow induces 
significant blockage and we see values of V' much less than I 
well beyond the tip radius. 

As seen above, some of the high speed tangential flow near the 
gear face plane is diverted into the tooth passage and a stagnation 
region appears, where this flow impinges on the leading surface, 
near the gear face. The net torque due to pressure effects can be 
represented by the difference between the leading and trailing 
surface pressure coefficients (C/.=(p-p eo )/0.5pV / - ip ), which is 
shown in Fig. 9. Clearly, the net torque is dominated by the near- 


DeloC. 



Fig. 9 A C P between gear leading and trailing tooth surfaces 
(one-half of symmetrical gear shown) for Diab gear 1 at 850 
rad/s 
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Fig. 10 Pressure torque per unit width versus axial coordinate 
(one-half of symmetrical gear shown) for Diab gear 1 at 850 
rad/s 



Machine Axis 


Fig. 11 Four notional shroud configurations for Diab gear 1 
geometry. Figure is to scale. The gray region defines the gear. 
The solid black lines indicate the positions of the large-axial 
and large-radial shroud walls. The dashed lines indicate the 
positions of the small-axial and small-radial shroud walls. 


face impingement observed in Fig. 6, In Fig. 10. the torque per 
unit width is plotted versus distance from the gear face. Here, 
torque is nondimonsionalized as 


n.vi- - 


f 




d\ 




( 7 ) 


where A/> is the pressure difference between the grid faces on the 
leading and trailing surfaces (which have identical i-r vertex co- 
ordinates). d.\ „ is the tangential projection of the area of the grid 
face, r is the radial coordinate of the grid face centroid, and ref- 
erence length L, cf is the gear tip radius. We see that, indeed, it is 
the near-face region that dominates the pressure windage loss 
torque. 


3.2 Gear Shrouding. As mentioned in the introduction, sev- 
eral researchers measured improvements in gear windage loss per- 
formance when shrouds of various configurations were employed 
[ 1-3.5]. In tins section, the aerodynamics of shrouding in die con- 
text of geometrically simple spur gear is studied and several con 
elusions are drawn, which may impact design considerations. The 
Diab gear 1 configuration is studied with four shrouding arrange- 
ments. nominally, large -axial-large -radial. small-axial-Iarge- 
radial, large-axial -small-radial, and small-axinl-small -radial Table 
2 and Fig. 1 1 quantify and illustrate the shroud dimensions, re- 
spectively. These four shroud configurations were chosen to be 
representative of the extrema of the full-shroud NASA Glenn tests 
detailed further below. The smaller axial and radial shroud clear- 
ances gave rise to mesh topology constraints that, in turn, led to 
poor quality block-structured meshes, unless overset meshes are 
used, as shown in Fig. 12. 

In Fig. 13. the predicted windage loss versus rotation rate is 
plotted for the unshrouded ease, validated and studied above, and 


Table 2 Diab shroud clearances 


Minimum 

Maximum 

Axial ( ,'ff U| 1 

0.0044 

0.1733 

Radial Off,, p ) 

0.0044 

0.0970 




Fig. 12 Cross section of overset mesh topology for the small- 
radial shroud 



Fig. 13 Comparison of predicted windage losses between free 
spinning Diab gear 1 and four shrouded configurations 
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Fig. 14 Comparison of predicted windage losses between four 
shrouded Diab gear 1 configurations 


the four shrouded configurations. All torn' of the shrouds give rise 
to very significant improvements in windage losses. The large- 
axial-large-radial shroud provides a 68% decrease in loss at 850 
rad/s: the small-axial-small-radial shroud provides approximately 
an 81% decrease at this speed. Indeed, distinguishing the perfor- 
mance gains between the shrouded cases is facilitated by exclud- 
ing the free spinning results, as in Fig. 14. Examining this plot, we 
see that reducing the axial and radial clearances from large to 
small provide approximately the same level of additional benefit 
over the large-axial-large-radial case, with the reduction of the 
radial clearance, providing somewhat more benefit in this particu- 
lar case. Applying both clearance reductions together provides the 
maximum benefit. 

Figure 15 shows the same pressure coefficient plot reported 
above but here, for the large-axial-large-radial shroud case, it is 
again at 850 rad/s. A highly edge loaded A Cp between leading and 
trailing surfaces can be observed as with the unshrouded case but 
with a very reduced range of A Cp. which, of course, leads to 
reduced torque. Figure 16 shows contours of pressure and surface 
shear stress lines on the leading surface for the four shrouded 
cases at 850 rad/s. It can be seen that the general features of 
impingement/stagnation near the face remain present for all of the 



Large Axial- 
Large Radial 


Large Axial- 
Small Radial 


Small Axial- 
Large Radial 



Small Axial- 
Small Radial 


Fig. 16 Contours of pressure and surface shear stress lines 
on the leading surface for the four shrouded cases at 850 rad/s 


configurations. 

The torque per unit width is plotted versus distance from the 
gear face for all four shrouds at 850 rad/s in Fig. 17. All four 
configurations exhibit edge loaded profiles, with their integral 
consistent with the net loss trend reported in Fig. 14. Comparison 
with Fig. 10 shows the dramatic reduction in pressure torque that 
has been achieved for these configurations but with retention of 
the basic features of the pressure torque distribution. 

An exploration of viscous losses was performed on the 
shrouded Diab gears. Figure 5 illustrated that the viscous losses 
per unit span for the unshrouded case increased with rotation 
speed and radius. The viscous losses per unit span are plotted for 
each of the four shrouded cases and the unshrouded case at 850 
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Fig. 15 A C P between gear leading and trailing tooth surfaces 
(one-half of symmetrical gear shown). Diab gear 1, large-axial- 
large-radial shroud, 850 rad/s. 



Fig. 17 Plot of torque per unit width from the gear face for the 
four shrouded cases 
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Fig. 18 Comparison of viscous torque per unit span for Diab 
gear 1 at 850 rad/s for the unshrouded and four shrouded 
configurations 


ratl/s in Fig. 18. It is observed that face shear is reduced for all of 
the shrouded cases due to the "couette-fiovvlikc" rotational bound- 
ary layer that arises between the gear face and the axial shroud. 
Also, of interest in this It cure are the greater viscous losses for the 
small axial shrouds compared with the large axial shrouds, rhese 
losses are still less than the unshrouded ease but suggest that 
viscous losses can increase as a percentage of total loss for very 
small axial shroud clearances. 

4 NASA Glenn Gear Windage Test Facility 

Hie NASA Glenn Research Center has recently installed a gear 
windage test facility [5]. The test facility is designed to param- 
etrize the effects of gear geometry, shroud geometries and sizes, 
different lubrication system configurations, system pressures and 
temperatures, and gear meshing cm windage loss. A sketch of the 
test facility is shown in Fig. 19. The facility has a 112 kW (150 
lip) dc drive motor that is connected to a 5.7:1 speed-up gearbox. 
The output of the speed-up gearbox is then connected to a torque 
meter prior to a coupling connection of the input shaft to the test 
gearbox. I lie input and output shafts have hydraulically operated 
clutches that allow the facility (in single or dual shaft mode) to be 
disconnected from the power applicr and/or magnetic brake at 
tadied to the output shaft. With the speed capability of the drive 



Fig. 19 Overview sketch of the NASA Glenn Research Center 
Gear Windage Test Facility [5] 
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Table 3 Basic gear dimensions [5] 


No. of irvih 

Module/pitch, mm (1/in.) 
Pace width. mm (m.J 
Pitch diameter, mm (in.) 
Pressure angle* deg 
Outside diameter, miu (in.) 


J2 

63 } (4) 
28.4 (1.12) 
330.2 (13.0) 
25.0 

342.65 (13.49) 

Table 4 Shroud wall clearances studied [5] 


Minimum 

Maximum 

Axial, mm (in.) 

0.762 (0.030) 

29.718 (1.17) 

Radial, mm (in.) 

0.762 (0 030) 

16.64 (0 655) 


motor and speed increasing gearbox and the dimensions of the test 
specimen, (he pitch line velocity can he taken to 280 m/s (55,00(1 
ft/min). 

The test gears can he run in the gearbox with or without shroud- 
ing. The shrouding clearance can be adjusted both radially and 
axially. This work reports data from lire first series of experiments: 
air-only, single-gear, and fully enclosed axial and radial shrouds. 
Details of tire seal- under test in this report are given in Table 3 
and the maximum and minimum shroud clearances are provided 
m Table 4 The gear has some modesl geometric complexities 
compared with tire idealized Diab gear studied above. These in- 
clude chamfered teeth, filleted teeth roots, and narrower body 
width between the teeth and the hub This work reports experi- 
mental data (’ll) results for the four extreme configurations 
large-axial-large-radial, large-axial-small-radial, small-axial large- 
radial. and small-axial-small-radial. A sketch of the shrouding ar- 
rangement is shown in Fig. 20. 

Data from the NASA Glenn Gear Windage Test Facility were 
measured in the following manner. Speed data arc measured using 
inductive pickups that read a (id-tooth disk on the end of each of 
the shafts. The output from the sensor (pulse/sine wave) is sent to 
a frequency to voltage converter. The output from the converter is 
then sent to a National Instruments card and read by l.AUVlbW. 
Data were taken at 10 Hz. The facility was operated at a series of 
increasing drive motor speeds. At each of these conditions, several 
different data were collected. The drive motor speed torque ap- 
plied to rotate the lesl hardware, internal shroud (fling-off) tem- 
perature. and internal shroud static pressure data were taken at 
steady drive motor speed conditions. Data were taken and then tile 
speed was incremented up a given drive motor speed from 500 



Fig. 20 Shroud assembly for test facility [5] 
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Large-Axial, Large-Radial 
Small-Axial, Large-Radial 
Large-Axial, Small-Radial 
Small-Axial.Small-Radial 


(Exp.) 

(Exp.) 

(Exp.) 

(Exp.) 



600 000 

Speed (radfs) 


1200 


Fig. 21 Experimental [5] and CFD results for the 13' pitch di- 
ameter NASA spur gear 


rpm to .’125 rpm (or 2587-16.168 rpm of the gem - shaft). The data 
were taken from motor controller speed, a commercially available 
torc|ue meter for torque, thermocouple for inside shroud tempera- 
ture (fling-off), and manometer. In order to determine the effects 
of the gear-only windage, a separate test was conducted with the 
entire system in place minus the gear. 

5 Experimental and Computational Results for the 
Glenn Facility Tests 

Grids were developed for the 13 in. (330.2 mm) pitch diameter 
spur gear configuration using similar topologies as the Diab cases. 
Four shroud configurations have been selected for study. Table 4 
lists the four baffle configurations studied computationally. 

Figure 21 shows comparisons of the NASA Glenn data and 
ovi-R-Ri-t predictions for the four shroud configurations. The data 
show Lbe same trends as the idealized Diab gear CFD studies 
reported above. Specifically, the large-axial-large-radial shrouding 
exhibits the highest loss levels, and the small-axial-small-radial 
shrouding exhibits the lowest loss levels. The benefit realized hy 
reducing both clearances is somewhat more substantial than for 
the Diab case. The CFD results are seen here to provide fairly 
good agreement with measured values. An interesting observation 
in the CFD results is that the large- axial -small-radial and small- 
axial-large-radial results are nearly identical along the entire speed 
line. 

The qualitative correspondence between the NASA Glenn and 
idealized Diab cases presented earlier suggests that die same 
physical loss mechanisms are acting. Figure 22 shows a view of 
predicted relative streamlines colored by pressure for the large- 
axial-large-radial shroud case at 700 rad/s. This image exhibits 
several of the features observed for the idealized case, including 
diversion of the near-face flow into the tooth passage and im- 
pingement on the leading surface, strong axial secondary vorticily 
in the tooth passage, strong ejection of this flow near the tooth 
centerline, and radial flow of the near-face streamlines below' the 
base of the teeth. 

6 Design Alternatives 

It was demonstrated above that axial and radial shrouding can 
reduce windage losses. Some ol the physics of these loss reduc 
lion schemes were studied there. Despite the experimentally and 
computationally observed differences in loss magnitudes between 



Fig. 22 Three-dimensional relative frame streamlines colored 
by static pressure for the NASA 13" pitch diameter spur gear, 
large-axial-large-radial shroud, at 700 rad/s 


unshronded and various shrouded configurations, in all cases, a 
significant component of the torques associated with windage 
arose from impingement onto die leading surface of the high ve- 
locity relative frame flow drawn into the tooth passage. Accord- 
ingly. in this section, we Teturo to the geometrically idealized Diab 
gear I in the large-axial -large-radial shroud configuration, and ex- 
periment numerically with four proposed tooth geometry modifi- 
cations aimed at mitigating this impingement and attendant wind- 
age torque. Figures 23—27 show an oblique view of the baseline 
geometry and the four alternative geometries considered: (1 ) lead- 
ing surface tooth-edge rounding. (2) leading+trailing surface 
tooth-edge rounding. (?) double slots on the top of the teeth, and 
(4) trailing surface ramp. 

Figure 28 shows a comparison of these lour simulations with 
the baselme large-uxial-large-radial case. The leading surface 
rounding and double slot geometries return nearly identical wind 
age loss. The leading -I- trailing surface rounding returns somewhat 
higher loss. However, the net loss obtained using the nailing sur- 
face ramp is approximately 30% lower than tire baseline configu- 
ration. The torque per unit width for the five geometries are plot- 
ted in Fig. 29. There, it can be seen that the ramp configuration 



Fig. 23 Predicted surface pressure coefficient and skin fric- 
tion lines for baseline tooth geometry 



Fig. 24 Predicted surface pressure coefficient and skin fric- 
tion lines for tooth geometry alternative: leading surface 
rounding 
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Fig. 25 Predicted surface pressure coefficient and skin fric- 
tion lines for tooth geometry alternative: leading+trailing sur- 
face rounding 


exhibits much smaller torques within the tooth channel and this 
clearly results in the reduced integrated loss for the entire gear. 

To further elucidate the physics involved in these numerical 
studies, we return to Figs. 24—27. In each of these figures, pre- 
dicted surface pressure coefficient contours are plotted along with 
selected surface skin friction lines. The baseline and two rounded 
geometries exhibit largely same qualitative flow features, with the 
rounded cases "smearing" the leading face impingement and trail- 
ing face detachment gradients. The tooth slots were conceived to 
“flush” the peak axial vorticity/low pressure regions with higher 
velocity incoming relative flow, thereby reducing the axial pres- 
sure gradient and thereby diverting the high relative velocity near- 
face flow into the passage. This appears to not have achieved the 
desired result. The trailing surface ramp geometry did have a sig- 
nificant impact on the aerodynamics. Specifically, the relative flow 
near the face is turned away from the gear. This turning induces a 
local pressure rise on the ramp, which contributes to windage 
torque. However, this flow has been diverted away from the tooth 
enough that subsequent diversion of this flow into the tooLh pas- 
sage has been virtually eliminated, resulting in almost no pressure 
rise on the leading surface. This gives rise to much less net torque 
shown in Fig. 29. and so the improved net performance of the 
ramp configuration observed in Fig. 28 is clearly due to the re- 



Fig. 26 Predicted surface pressure coefficient and skin fric- 
tion lines for tooth geometry alternative: tooth slot 



Fig. 27 Predicted surface pressure coefficient and skin fric- 
tion lines for tooth geometry alternative: trailing surface ramp 

Journal of Fluids Engineering 



Fig. 28 Windage loss predictions for the baseline tooth geom- 
etry (Diab gear 1 with large-axial-large-radial shroud) and four 
geometric alternatives 


duced integrated tooth passage torque more than offsetting the 
increased torque associated with the ramp turning itself. 

7 Conclusion 

This paper has summarized a number of CFD studies focused 
on the aerodynamics of gear windage losses, lire goals of this 
work have been to validate the numerical and modeling ap- 
proaches used for these applications and to develop physical un- 
derstanding of the aerodynamics of gear windage loss. Compari- 
sons are made with experimental data from open literature, and 
data recently obtained in the NASA Glenn Research Center Gear 
Windage Test Facility. Parametric shroud configuration studies 
carried out in the Glenn experiments and the CFD analyses eluci- 
date important physical mechanisms of windage losses as well as 
mitigation strategies due to shrouding and newly proposed tooth 
contour modifications. 

The following conclusions apply. (1) The NASA Glenn data 
and CFD analyses show that axial and radial gear shrouding are 
effective in significantly reducing gear windage losses both inde- 
pendently and when employed together. (2) For all configurations 
studied, the dominant physical mechanism contributing to wind- 





Fig. 29 Torque per unit width predictions for the baseline 
tooth geometry (Diab gear 1 with large-axial-large-radial 
shroud) and four geometric alternatives at 850 rad/s 
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age losses is the pressure field associated with diversion and im- 
pingement of the high speed relative flow on the leading tooth 
surface. (3) Shrouding mitigates the magnitude of this mechanism 
but not its dominance in the loss budget. (4) Although viscous 
effects have been identified as a secondary source of windage 
losses, they are not insignificant. Also, it was identified here and 
in Refs. [15,16] that the sublayer resolved turbulence model used 
here is still more accurate than a high Reynolds number form with 
wall-functions but apparently underpredicts viscous losses for the 
rotational boundary layer studied. Accordingly, future attention to 
turbulence modeling is important. (5) Viscous losses are reduced 
by shrouding as well. However, small axial clearances were ob- 
served computationally to increase viscous losses compared with 
larger shroud clearances, suggesting that this effect could become 
important for even smaller clearances and at higher speeds. (6) 
The CFD results show good agreement with the Glenn experi- 
ments. (7) The NASA Glenn data show quite similar loss trends to 
the idealized shrouded gear configurations studied computation- 
ally. (8) The CFD studies suggested a set of possible geometric 
tooth modifications to further reduce windage loss. One of these 
appears quite promising and the authors anticipated studying this 
in the NASA Glenn facility. 

Our continuing work in the gear windage aerodynamics area 
focuses on (1 ) further validation for shrouded gears as more Glenn 
data become available, (2) validate for viscous-only losses with 
available data on rotating disks [4,23], (3) multiphase flows, (4) 
overset gridding in the context of gear tooth contact, and (5) fur- 
ther evolution of design guidance for minimizing windage losses. 
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volume 
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Ill the final stages of the NRA program, we turned our attention to further physics 
interpretation/understanding, and the role of multiphase flow. We discovered and proved that angular 
momentum confinement near the pitch-line-face region plays an important role in shroud effectiveness. 
This was a discovery not understood by the community prior to our research. As postulated by other 
researchers, we showed that tooth channel mass flow rate restriction due to the shrouds is also important 
in reducing windage loss, and is the dominant mechanism for very small clearance radial shrouds. It was 
demonstrated that even relatively small leakage paths in the shroud can reduce the shroud’s 
effectiveness due to the loss of angular momentum confinement. 

In the multiphase flow modeling area, a plausible set of droplet and film models, separately 
validated for aerodynamic application, returned significantly increased losses associated with oil 
lubrication. The multi-phase studies undertaken demonstrated a richness of physics that we have only 
begun to explore, including significantly larger predicted losses than suggested by homogeneous 
assumptions, important local droplet concentration effects, and the necessity of incorporating film 
modeling. The 2012 AHS Forum paper which summarizes this research is included here in its entirety. 
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ABSTRACT 

This work, builds on the authors’ recent CFD research on high speed gearbox windage loss^ 1 ' 4 '. There, 3D steady and 
unsteady, single phase simulations of spur gears were presented. Those studies served to elucidate some of the physics of 
windage losses, validate the methods and models employed against available experimental data, and explore design 
alternatives to mitigate windage loss. With this validated capability in hand, we have recently studied computationally two 
critical design relevant components of high speed gearbox aerodynamics: 1) The role of shrouds, widely used in rotorcraft 
transmission systems for windage loss mitigation, and 2) the role of lubricant transport, i.e., the multiphase flow resident in 
real gearboxes. Research findings and attendant recommended design guidance are reported in this paper. 
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inter- field non-drag force kernel 
mass flow rate 
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mean kinetic energy 
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Weber number 
volume traction 
internal energy, turbulence dissipation rate 
molecular viscosity 
density 
surface tension 
inter-field mass transfer rate 
system rotation rate 
viscous stress tensor 


V 


control volume 


INTRODUCTION 

Gear windage losses refer to the power losses that are a 
result of aerodynamic forces (pressure and viscous) acting 
against the rotation of high speed gears. Windage losses are 
a source of significant heating and power consumption in 
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rotorcraft and other gearbox systems where high pitch-line 
velocities are encountered (nominally >80 m/s). 

In earlier publications by the present authors 1-4 , we 
studied various aspects of windage losses, using CFD, in the 
context of single phase flow, and validated our method 
against open literature data 5 , and data recently obtained in 
the Gearbox Windage Facility at NASA Glenn Research 
Center 6 . That work showed that the dominant physical 
mechanism contributing to windage losses in spur gears is 
the pressure field associated with diversion and impingement 
of the high speed relative flow on the leading tooth surface 
adjacent to the face. Shrouding the gear was observed to 
mitigate the magnitude of this mechanism but not its 
dominance in the loss budget. Also, this mechanism 
suggested some passive gear geometry treatments to further 
reduce windage loss, these were explored computationally 
and one concept has shown design promise 3 . This earlier 
work complements a number of CFD studies carried out by 
other researchers that have appeared in recent years 7,8,9,10 , 
some of which are reviewed in Ref. 3. 

In this paper we report on progress in two research 
tracks. Firstly, a computational exploration of the physical 
mechanisms associated with shrouding has been carried out 
and several conclusions are drawn related to “why shrouds 
work” and “when/why they do not”. Secondly, we 
incorporate non-equilibrium two-phase flow analysis and 
show that droplet and film flow distributions play an 
important role in windage loss. 

It was shown experimentally by Dawson 11,12 and 
Winfree 13 that geometric modifications of the near-gear flow 
path, such as shrouding, can significantly reduce both 
windage losses and lubricating oil consumption (80% and 
40% reductions observed, respectively, in Ref. 13). 
Accordingly, shrouding is of significant interest in rotorcraft 
transmission design and is used in numerous systems today. 
Unfortunately, only a few empirical studies have been made 
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available to date, and little physical understanding has been 
proposed. Indeed, Winfree 13 observed that shrouding a bevel 
gear reduced the ability of the gear to pump fluid from the 
gear axis and through its teeth, which suggests that the 
shroud acts to restrict the amount of fluid that the gear draws 
in and then expels through the tooth channel. Dawson 11 ' 12 
also attributed windage losses in spur gears to the drawing of 
air axially into the tooth channel and its ejection radially at 
the gear centerline, i.e. that the gear acts as a pump and that 
windage losses can be interpreted as the rate of work 
performed on the air in this pumping process. He supported 
this hypothesis with flow visualization, and an observed 
reduction in windage losses when the channel entry was 
blocked - as by an axial or radial shroud. These hypotheses 
are supported in part below, but the mass flow blockage 
mechanism itself is shown to be only part of the story, with 
angular momentum dynamics playing importantly as well. 

Accordingly, there is a need for more experimental data 
and analysis aimed at improved understanding of the 
physical mechanisms involved in shrouding schemes for loss 
reduction. The Gearbox Windage Facility at NASA Glenn 
Research Center 6 , has been established in part to produce 
such data, and we recently validated the CFD method used 
here against a number of shroud configurations studied in 
that facility 3 . Below we study several shroud configurations 
and formulate and support a hypothesis as to why they are 
effective. These studies also show some instances where 
they are not effective. Together these observations lead to 
design guidance proposals. 

The important role of lubrication oil in windage losses 
is well known although not well understood. Most high 
speed experimental windage studies to date have been 
performed in air 5 ' 6 ’ 13 , though forthcoming results from the 
NASA Glenn Windage Facility will provide 2-phase 
validation data. The incorporation of the two-phase gearbox 
environment in predictive analysis is generally quite 
primitive; all of the correlations in the literature employing, 
at most, dimensionally consistent mixture density and 
viscosity parameters 8 , obtained using estimated mass loading 
and droplet homogeneity assumptions. Imai et al. 10 are the 
only group to the authors’ knowledge to perform multiphase 
CFD analysis for gearbox aerodynamics. They analyzed a 
two-gear spiral-bevel system using a homogeneous VOF 
method and porosity model for the two-phase flow and gear 
meshing, respectively. These authors demonstrated that the 
VOF method was able to capture details of the macro- 
physics of lubrication jet breakup and transport, but such a 
method is inherently incapable of capturing the details of 
droplet/mist formation and transport at the scales where they 
contribute to mass loading of the environment (l-300pm). 

In our view, the minimum level of modeling required to 
capture two-phase environmental effects in high speed 
gearboxes are methods that account for non-equilibrium 
interfacial dynamics. That is, two-fluid Eulerian or 
Lagrangian methods can accommodate slip velocity between 
the droplets and carrier air through the modeled influence of 


(at a minimum) interfacial drag. The reason that a non- 
equilibrium approach is needed is that droplet distribution is 
fundamental to local mass loading which in turn defines 
local mixture density and viscosity, and thereby loss 
magnitudes directly. The large shears and centrifugal forces 
encountered in high speed gearboxes lead to droplet 
distributions that are highly non-uniform by virtue of the 
droplet size spectrum that arises (i.e. the majority of droplet 
mass would have to reside in droplets less that 1pm in 
diameter for a near homogeneous distribution to arise). 
Lagrangian methods can capture these physics, but are likely 
to be very computationally intensive for the gearbox 
problem, and inherently unable to handle the wall films that 
arise. Accordingly, we have chosen an Eulerian two-fluid 
approach, which can be instrumented to model all of the 
physical mechanisms present in these systems. Specifically, 
it enables the modeling of droplet dynamics (drag, 
dispersion), droplet wall impaction and film breakup (both 
inter-field mass transfer mechanisms), multiple droplet sizes, 
droplet breakup and coalescence, film shear, and potentially 
other mechanisms. Most of these are implemented below. 

Below we study a shrouded spur gear configuration 
under multiphase flow conditions. Oil mass loading, droplet 
size, and number of droplet fields are parameterized. Also, a 
scheme to properly accommodate wall films is introduced. A 
number of conclusions are drawn from these studies related 
to mass loading specification, droplet distributions, the 
importance of modeling films and boundary conditions. In 
particular it is shown that simple mixture density and 
viscosity based loss correlations are unlikely to return 
reliable quantitative results. 

This paper is organized as follows. First, the governing 
equations, models and numerical procedures are reported. A 
series of shrouding related numerical studies are then 
performed and discussed. This is followed by a series of 
multiphase flow analyses, which demonstrate some of the 
richness of the multiphase physics that arise in these 
systems, suggest shortcomings likely to arise in simpler 
analysis, and point to further modeling refinements to be 
considered. 

THEORETICAL FORMULATION 
Governing Equations 

A two-fluid Eulerian formulation is adopted. The 
conservation of mass, momentum and energy for each 
constituent present are written in integral conservation law 
form for a compressible flow in a steadily rotating 
coordinate system, in Equations 1-3. k, 1 are constituent 
indicators, here k — 0, 1, 2 for air, droplets, and films 
respectively. For all simulations presented in this paper, the 
carrier air field is assumed a perfect gas; the disperse droplet 
and continuous film oil fields are treated as mcompressible. 

J p k a k v .dS = Y,(y‘ k -T u ) (1) 
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Turbulence and Multiphase Modeling 

A high-Reynolds number k-s turbulence model with 
wall functions is used in the studies that follow, for the air 
field only. Although sublayer resolved turbulence modeling 
can improve the accuracy of the viscous component of loss 
in gear windage analysis 1,3 , its underprediction by high 
Reynolds number formulations are of secondary importance 
in the shroud and multiphase flow studies undertaken here. 
Multiphase turbulence effects are discussed below. 

The terms D a , I M and A fi 1 that appear in the continuity 
and momentum equations represent interfacial drag forces, 
interfacial mass transfer and interfacial “non-drag” forces, 
respectively. A broader discussion of the relevance of these 
terms in general two-fluid simulations, as well as details of 
their numerical implementation is available in Ref. 14. Here 
we present the models used for the gearbox windage 
application. 


There is no data yet available with which to calibrate 
droplet size distribution in the gearbox system. In this first 
multiphase flow study, we account for droplet size by 
parameterizing its effect on the solution. Droplet 
coalescence and breakup are not yet modeled. 

In the high speed gearbox system, a significant (but 
unknown) amount of the oil in the gearbox is constituted as 
films on the gears, shrouds, housing and other components. 
Although 2-field analysis (air+droplets) can provide some 
interesting physics and design guidance for relevant mass 
loadings when shrouds are present, droplets will accumulate 
locally in a non-physical manner if not allowed to deposit as 
films on the walls since their interfacial area density and 
thereby drag is so high. Modeling the mass transfer from 
droplets to films enables films to be modeled explicitly with 
relevant (small) interfacial areas densities, and appropriate 
drag kernels. For droplet deposition to films we incorporate 
a hybrid impaction and diffusion mass transfer model due to 
Haworth et al. 17 : 
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For drag exerted on the (assumed thin) films by the air 
we have: 


Droplet drag is modeled in a conventional manner 1 " : 
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To accommodate the effect of the carrier field (air) 
turbulence, a droplet dispersion model due to Lopez de 
Bertodano 16 is used: 


= C TD p a,r k s,r V a drop, “ (5) 

The drag and dispersion models in Equations 4 and 5 
represent a “minimum” set of interfacial force models in that 
for low mass loading, they are sufficient to account for the 
distribution of droplets within the gearbox. We do not model 
other interfacial dynamics in this work (e.g., lift, virtual 
mass) as they are likely to be negligible in this high density 
ratio system (p ai 7p drc,ple ‘ « 1). 
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Lastly, we model film breakup mass transfer to droplets 
(in order that a physically representative droplet-film 
equilibrium is attained) following Schadel and Hanratty 18 : 
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Equation 9 is applied when the local Weber number exceeds 
tlie critical breakup Weber number of 17. 


CFD Numerics and Code 

The CFD code used in this work, NPHASE-PSU h is a 
parallel face-based, cell-centered, arbitrary-element 
unstructured multiphase flow solver. The baseline algorithm 
follows established segregated pressure based methodology. 
Compressibility is accommodated by including a density 
correction in the formulation of the pressure corrector 
equation. A number of algorithmic details associated with 
multiphase flow are relevant here, and those details are 
available in Ref. 14. The code is instrumented with overset 
meshing and mesh motion, both used in our earlier windage 
publications, but all of the simulations carried out here use 
contiguous meshes which are stationary in the rotating frame 
of reference. A single tooth passage is simulated and 
periodic boundary conditions are applied. Since these are 
spur gears the symmetry of the system was also exploited 
and only 14 of the domain was modeled. A 517,348 cell 
mesh was used for the unshrouded analyses presented below, 
fewer cells are required for the various shrouded cases. The 
meshes conform to wall spacing requirements for the 
turbulence model wall functions used (y + <200). These grids 
are consistent with those used in our earlier work 1 '. 


All of the studies presented in this paper are applied to 
the same spur gear geometry, a 72-tooth spur gear studied 
experimentally by Diab et al. 5 , which has a pitch diameter of 
288mm (r tt p=1.48m), a tooth width of 30 mm and a module 
of 4 mm. 


SHROUDING STUDIES 

It has been known for some time that axial and radial 
shrouds can significantly reduce gearbox windage 
losses 11 ' 12,13 . In Refs. 2-4 it was shown that the favorable 
impact that shrouding has on windage losses can be 
predicted quantitatively, including the trends associated with 
the axial and radial clearance of the shrouding. However, the 
physics behind why shrouding schemes work has remained 
only partially understood. 

It has been conjectured in the literature 11 ' 12,13 that 
shrouds restrict the “pumping” ability of the gear teeth, and 
that this is fundamentally responsible for the reduction in 
windage loss. To explore this, we study unshrouded and 
shrouded spur gears computationally. Figure 1 shows a 
visualization of the unshrouded simulation of the test spur 
gear rotating at 1000 s 1 in the direction indicated by the 
arrow. There, relative frame streamlines show' the drawing in 
of the near face flow and radial ejection of this flow at the 


tooth channel centerline, consistent with Dawson’s 
observations 11,12 . A computational shce is taken across the 
gear top land and contour s of radial velocity are plotted with 
the ejecting relative flow streamlines. When integrated, this 
net radial mass flow rate provides a quantitative figure for 
the gear “pumping”. 

The four shrouds studied are designated large -radial- 
large-axial (LRLA), large-radial-small -axial (LRSA), small- 
radial-large-axial (SRLA), and small-radial-small-axial 
(SRSA). These are shown schematically in Figure 2. The 
small and large radial shrouds extend .0044 ’T tip and 
. 1733*4(1,, radially outward from the tip. The small and large 
axial shrouds extend .0044 *r„p and ,0970*r tlp axially away 
from the gear face. 



Figure 1. Relative streamlines and radial velocity 
contours for an unshrouded simulation of a spur gear. 



Figure 2. Predictions of net mass flow ejected radially 
from gear tooth channel for unshrouded and four 
shrouded configurations. 
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Figure 3. Windage loss predictions for shrouded and 
unshrouded spur gears. 

The radial mass flow integration was carried out for 
four shroud configurations and the results plotted in Figure 
2. Figure 3 shows predicted windage loss vs. rotation rate for 
the unshrouded case and the four shrouded cases. The 
unshrouded simulation compares well with the experimental 
measurements due to Diab et al. 5 , and our previously 
published results using another CFD code 3 . (The present 
results somewhat overpredict the loss compared to 
experiment and the earlier analyses - the reasons for this 
have not been studied, but these differences do not impact 
the conclusions drawn from the comparative studies 
undertaken. ) 

Comparing Figures 2 and 3, it is clear that the integrated 
mass flow rates do not exhibit the same relative trends as the 
loss. Specifically, considering the LRLA shroud, the 
reduction in mass flow from the unshrouded gear, =25%, is 
far less than the observed =72-74% reduction in windage 
loss. For the SRSA shroud, the reduction in mass flow from 
the unshrouded gear, =97%, is far more than the observed 
=90% reduction in windage loss. In Figure 4, the ratio of 
windage loss to integrated mass flow is plotted for all cases. 
From this figure we see that each configuration exhibits 
similar trends with gear speed, but that the different shroud 
schemes exhibit large variations in this ratio. Indeed, the 
LRLA and SRSA shrouds show a consistent order of 
magnitude difference in this parameter. So clearly, mass 
flow restriction is only part of the physics which constitute 
performance gains with shrouds. 

In order to better understand the nature of the 
shrouding’s effect on losses, consider Figure 5a which 
shows contours of absolute tangential momentum for an 
unshrouded and the LRLA shrouded analyses of the test gear 
at a rotation rate of 1000 s' 1 . The contour slice is taken at 
mid channel. In all gear systems the gear acts to impart 
significant angular momentum to the fluid. It can be seen 


that the shrouded gear confines the region of large absolute 
angular momentum near the teeth whereas this momentum is 
transported away from the teeth in the unshrouded case. 



Figure 4. Predictions of loss/mass flow ejected radially 

from gear tooth channel, for unshrouded and four 
shrouded configurations. 

Figure 5b shows a plot of angular momentum along an 
axial line through the domain at the pitch radius (indicated in 
Figure 5a) for the unshrouded and all of the shrouded cases. 
Considering the LRLA case, it is observed that the 
difference in angular velocity between the teeth and the 
near-gear-face flow is much lower for the shrouded case due 
to the confinement of the angular momentum. This leads to a 
much lower deceleration of the flow that impinges on the 
leading tooth near the face, in turn leading to lower pressure 
in the stagnation region where the flow is drawn in, as 
observed in Figure 5c. So for the LRLA, we have a small 
effect on loss reduction due to mass flow restriction, and a 
larger effect due to angular momentum confinement; 
together these give rise to a net loss reduction across the 
speed range studied of between 72% and 74%. This budget 
of mass flow and momentum effects is reflected in Figure 4 
where the LRLA data is a factor of approximately 3 lower 
than the unshrouded budget and a full order of magnitude 
less than the SRSA budget. 

Consider next the SRSA case. Figure 5b shows that the 
angular velocity between the teeth and the near-gear-face 
flow is again higher than for the unshrouded case but is 
lower that the large-axial cases by virtue of the proximity of 
the axial shroud. So we expect a reduced benefit due to 
angular momentum confinement for the small-axial shrouds 
compared to the large axial shrouds. Also, it was observed in 
Refs. 2, 3, 4 that increases in viscous loss begin to manifest 
for small axial shrouds. So we have the situation where mass 
flow rate is dramatically reduced by some 97% for all 
speeds, but the losses due to reduced angular momentum and 
viscous effects lead to higher relative losses than the reduced 
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a) 


■Iimdi . 


V e (m/s) 

1 851e+02 
1 388e+02 
9.255e+01 
4.628e+01 
0.000e+00 



DeltaP (Pa) 

4 166e+Q3 
2 502e+03 
8.372e+02 
-8.273e+02 
-Z492e+03 


mass flow rate alone would suggest. This budget of mass 
flow and momentum effects is reflected in Figure 4 where 
the SRSA data is a factor of approximately 3 higher than the 
unshrouded budget and a full order of magnitude higher than 
the LRLA budget. 

Further consideration of Figures 4 and 5b lead to a final 
observation. The LRSA and SRSA shrouds exhibit similar 
near face angular momentum profiles as do the LRLA and 
SRLA shrouds. However, the SRSA shroud shows predicted 
windage losses 53-55% lower than the LRSA shroud across 
the speed range, and the SRLA shroud shows predicted 
windage losses 55-56% lower than the LRLA shroud across 
the speed range. This tracks well with the mass flow rate 
reductions: SRSA shows predicted mass flow rate reductions 
87- 88% lower than the LRSA, SRLA shows predicted mass 
flow rate reductions 83- 85% lower than the LRLA. So for a 
given axial shroud location we see that radial shrouding’s 
dominant effect is indeed on mass flow reduction. 

The roles of angular momentum confinement and mass 
flow restriction through the channel can provide significant 
guidance in the design of gear+shroud systems. In Refs. 2, 3, 
4, several gear tooth modifications were explored in order to 
reduce the mass flow into the channel, a spoiler concept 
shows promise there. Analogously, CFD can be used to 
optimize shroud geometry towards reducing channel mass 
flow anchor to maximizing angular momentum confinement. 
Considering the latter, in real gearbox systems, holes are 
required in the shrouding for plumbing, drainage, meshing 
and hub penetration. Such penetrations can reduce the 
effectiveness of shrouds by allowing angular momentum to 
“leak” away from the near-pitchline-face region. To explore 
this. Figure 6 shows elements of a drainage slot 
parameterization study. There, the results of three 
simulations are presented, with the nominal, fully enclosed 
SRLA shroud explored above, and two slotted-shroud 
configurations, each rotating at 750 s' 1 . 



Figure 5. a) Contours of absolute angular velocity for 
unshrouded and shrouded (LRLA) simulations of test 
gear at 1000 s' 1 , b) Angular velocity along pitch line for 
all configurations at 1000 s , c) Tooth pressure contours 
for unshrouded and shrouded (LRLA) simulations at 


1000 s 


Figure 6. a) Contours of absolute angular velocity for 
SRLA configuration at 750 s' 1 , a) Fully shrouded, b) 
Radial shroud corner slots, c) Radial shroud centerline 
slot 
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Slots were introduced at the radial shroud corners and 
centerline as shown in the figure, and the computational 
domain was extended radially and axially into a volume so 
that the flow can exhaust. Figure 6 shows predicted contours 
of angular momentum. Both slotted geometries allow some 
transport of angular momentum into the external volume. As 
seen in the contours, this comes at the expense of angular 
momentum in the near-pitchlme-face region, for both slotted 
cases, when compared to the fully enclosed run. Also, the 
ejected channel mass flow rate increased by 13% for the 
comer slot, and decreased by 10% for the center slot 
(perhaps counter-intuitively). Figure 7 shows the loss values 
associated with these simulations. For these notional slots, 
windage loss is approximately doubled. This significant loss 
increase associated with comparatively small shroud 
penetrations is consistent with data recently published by 
Handschuh and Hurrell 6 These simple parametric slot 
studies indicate that unavoidable penetrations in gear 
shrouds can significantly hurt the windage loss performance 
of these systems, and therefore CFD optimization of the size 
and placement of these penetrations can benefit gearbox 
performance. Of course, due to the non-periodic nature of 
real penetrations, such optimization studies will necessarily 
be much more computationally intensive than the 
exploratory studies presented here. 



Figure 7. Windage loss predictions for shrouded gears. 

Loss values for SRLA slotted configurations are plotted 
at 750 s' 1 . 

MULTIPHASE FLOW STUDIES 

Windage loss in lubricated systems is mtuitively 
dependent on droplet concentrations that are nominally local 
to the gear pitchline. Therefore CFD methods that model 
droplet concentration and size distributions accurately will 
be required in order that quantitatively accurate loss 
predictions are returned. This is a significant challenge for a 
number of reasons including: 1) There is no experimental 
data yet available with which to calibrate/validate droplet 


concentration or size distributions in the gearbox system, 2) 
Accurate and validated models for droplet drag, dispersion, 
and deposition, film drag and breakup are required - a first 
set of these models is incorporated here. But also, in order 
that local droplet size distributions be predicted accurately, 
models for droplet coalescence and breakup, and droplet 
bouncing and splash will need to be incorporated, as likely 
will more sophisticated models for the baseline set of 
mterfacial transfer models used here (Equations 4-9), 3) The 
boundary conditions for these multiphase simulations are 
difficult to define unless the complete gearbox with all lube 
injection plumbing, shroud features and gear meshing are 
modeled. Indeed, as discussed below, even estimating the 
total oil mass loading for the gearbox is challenging. 

In these first multiphase flow studies, we account for 
droplet size by parameterizing its effect on the solution, and 
by performing multiple-droplet- field analysis. The droplet 
models presented in Equations 4-6 (drag, dispersion, 
deposition) have been partially validated against a high 
speed aerodynamic data set 19 . Figure 8 shows predicted non- 
dimensional deposition rate vs. distance along the airfoil 
surface for an MS-317 airfoil at 0° angle-of-attack, free 
stream velocity of 90 m/s, and a free stream median 
volumetric diameter of 92 pm. In the experiment, 10 droplet 
size “bins” were used to characterize the far-field conditions 
- free stream liquid water content and droplet diameter were 
reported for all ten. Accordingly, a 12 field simulation was 
run (air + film + 10 droplet fields). Figure S shows that the 
NPHASE-PSU results match closely with the models in 
LEWICE-3D 20 , for the case where a droplet splashing model 
was not-mcluded. LEWICE does a much better job of 
matching the data when a splashing model is included, 
suggesting we include splashing in our future work. 


1 C LEWICE Splashing 

LEWICE No Splashing 

09- - NPHASEPSU No Splashing 



Figure 8. Non-dimensional deposition rate vs. distance 
along airfoil surface for an MS-317 airfoil. Contours of 
predicted total droplet volume fraction are also shown. 
Data and LEWICE predictions from Refs. 19, 20. 
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The single phase simulations presented above included 
no inflow or outflow boundaries. This is accommodated 
numerically by fixing the static pressure near atmospheric at 
a single cell in the computational domain, an approach that 
is guaranteed to conserve mass. In the multiphase analyses, 
this cannot be done since there is no mechanism to replenish 
mass loss errors in a given field during the course of 
iteration. Accordingly we introduce notional inflow and 
outflow boundaiy-’ conditions to the test gear configuration, 
and abandon the model symmetry. Specifically a small 
(compared to the pitchline velocity) inlet flow of air and 
droplets is specified and an outflow boundary is also 
specified, as shown in Figure 9. 



Figure 9. Computational domain for multiphase 

analyses. 

The first set of calculations performed employ one 
droplet field and no film field. Per Equation 1, we adopt a 
volume fraction form for constituent transport. The air and 
droplet volume fractions are specified at the inlet in addition 
to their velocities, and this defines the oil injection mass 
flow rate. However, although the injection flow rate 
influences the oil mass loading (= [total mass of oil]/[total 
mass of oil+air]) in the gearbox (or shr ouded gear system), it 
does not define it. Specifically, the mass of oil in the 
computational domain depends on many factors, but is not 
uniquely defined by inflow and outflow mass flow rate since 
the steady state oil continuity equation requires: 

J p°' l a 0, ‘ u* 1 -dA - 0 , (10) 

A 

but says nothing about the distribution of oil in the domain. 
Tins lias two implications: 1) Measuring lubn cation flow 
rates experimentally is not sufficient for validation - we will 
need to know local mass fractions inside the enclosure, 
2) Computationally we cannot “specify” mass loading, it is 
an output of fire computation. So although increasing inflow 


flux of oil to the rig (or computational domain) will increase 
the equilibrium oil mass loading in the system, the 
relationship will not be linear. 

For the 2-field droplet studies we ran a range of fixed 
droplet sizes, inlet oil flow rates, and gear speeds, all using 
the test gear with radial shroud extending .101*r bp radially 
outward from the tip and axial shroud extending .312*r tlp 
axially away from the gear face. Figure 10 shows the results 
of a cross section of these runs. Several observations are 
relevant: 1) All multiphase analyses exhibit a noticeable 
increase in windage loss over the single phase case, 
2) Increased mass loading for a given droplet size increases 
loss, 3) Increasing droplet size for a given mass loading 
increases loss. 

Pitch Line Velocity (m s ') 



Figure 10. Predicted windage losses for single phase and 
two-field air -droplet simulations. 

Here we interrogate the details of the simulations in 
order to better understand the physics underlying these 
windage loss trends, and to ascertain how general these 
results might be. In Figure 11a, contours of droplet volume 
fraction are plotted for a 2pm droplet field at a rotation rate 
of 750 s' 1 . The contour slice is taken at mid channel and the 
tooth region is displayed as partially transparent. Figure 11b 
shows the same plot for the 16pm droplet field. 

As expected, the 2pm droplet distribution is more 
uniform than the 16 pm case. Nevertheless, even for the 
smaller droplets, with their attendant high drag coefficient, 
the high shear rates, centrifugal forces and secondary air 
flows in this system give rise to a significant accumulation 
of droplets near the shroud and “starvation” in the tooth 
channel region. This effect is much more pronounced for the 
16 pm droplets which exhibit almost no droplets in the tooth 
channel region and even more accumulation of droplets near 
the shroud comers. In both cases a slight asymmetry 
associated with the non-symmetric boundary conditions is 
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noticeable. The predictions in Figure 10, 11a, and lib 
suggest that a simple mixture density scaling will not be 
suitable for loss correlations. Droplet diameter should be 
incoiporated in some fashion. Figures 12a and b isolate the 
low and high mass loading conditions respectively. Included 
are the single phase results for this shroud and the single 
phase result scaled by the bulk mixture density based on the 
mass-loading. Clearly in both cases the predicted losses are 
much higher than would be returned by simply scaling a 1- 
phase analysis by the bulk mixture density or, equivalently, 
by running a homogeneous mixture model. As expected, the 
2pm results are closer to the scaled 1 -phase loss than the 
16pm results, but they are still significantly higher, 
especially at high speeds. 


a) 



alpha_1 

1 000e-02 
1 000e-04 
1 000e-06 
1 .000e-08 
1 000e-10 



Figure 11. FTedicted droplet volume fraction distribution 
a) 2p.m. b) 16 (im simulations at 750 s' 1 . 



Figure 12. Predicted windage losses for single phase and 
two-field air-droplet simulations. Mixture density scaled 
1-phase results included, a) Low mass loading, b) High 
mass loading. 

Apparently, the non-equilibrium interfacial dynamics 
have led to non-uniformities in the mixture that increase 
windage losses. In order to better understand the physical 
reasons for this, consider Figure 13, where the higher mass 
loading, 750 s' 1 , 2pm and 16pm simulations are compared. 
Contours of mixture density (p m =Pair«air+P uictmi) are plotted 
on the gear symmetry plane. In the 2pm simulation, droplets 
accumulate near the advancing face, not only at the face 
edge, but near the tip all along the channel, leading to higher 
local mixture densities and pressures there, and thereby 
higher net Ap across the channel. The net predicted windage 
loss for this case was 515 W vs. 392 W for single phase. 
The 16pm case shows significant mixture density increase 
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all along the leading face, and a larger pressure magnitude 
than the single phase case, albeit qualitatively localized like 
the single phase case. The net predicted windage loss for this 
case was 1 165 W. 



Figure 13. Elements of higher mass loading, 750 s 1 
simulations, a) 2|tm case. Contours of mixture density at 
gear centerline (left) and tooth static pressure contours 
(right), a) 16jim case. Contours of mixture density at 
gear centerline (left) and tooth static pressure contours 
(right). 

As discussed above and evidenced in Figure 11, 
droplets will accumulate locally in anon-physical manner if 
not allowed to deposit as films on the walls since their 
mteifacial area density and thereby drag is so high. 
Accordingly we have implemented a film field which will 
exchange mass with the droplet field (droplet deposition, 
film breakup [Equations 6 and 9]) and has an associated air- 
film drag law (Equation 7) which accommodates the much 
lower interfacial areas density of the film and thereby allows 
it to flow along walls. The CFD code is instrumented with 
logic to allow films to be resident only in wall adjacent cells. 
We present preliminary results in Figure 14, showing the 
film volume fraction for a 3-field, 2|±m, 750 s' 1 , low mass 
loading run. Films develop along the walls, with maximum 
volume fraction (i.e., thickness) on the leading face and 
shroud corner, and mimima near the channel centerline. 


Droplet volume fraction contours for the same case shows 
that the droplets exhibit maximum concentrations near the 
shroud comers, but do not accumulate locally in a non- 
physical manner (compare Figure 11), since they are locally 
deposited as film which is transported, thickens and 
exchanges mass again with the droplet field along the axial 
shroud. These admittedly preliminary results exhibit the 
desired/required droplet film behavior as we move forward 
with model refinement and validation. 



f 


Figure 14. Elements of 3-field, 2|im droplet, lower mass 
loading, 750 s 1 simulation. Contour- s of a) Film volume 
fraction, b) Droplet volume fraction. 
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CONCLUSIONS 

This paper has summarized a number of CFD studies 
focused on the aerodynamics of high speed gear windage 
losses. The goals of this work have been to explore the roles 
of shrouds in these systems towards improved physical 
understanding and design guidance. The following 
conclusions apply: (1) Angular momentum confinement near 
the pitch-line-face region plays an important role in shroud 
effectiveness. (2) As postulated by other researchers, tooth 
channel mass flow rate restriction due to the shrouds is also 
important in reducing windage loss, and is the dominant 
mechanism for very small clearance radial shrouds. (3) Even 
relatively small leakage paths in the shroud can reduce the 
shroud’s effectiveness due to the loss of angular momentum 
confinement. (4) A plausible set of droplet models, 
separately validated for aerodynamic application, return 
significantly increased losses associated with oil lubrication. 
(5) The preliminary multi-phase studies undertaken 
demonstrate a richness of physics that we have only begun 
to explore, including significantly larger predicted losses 
than suggested by homogeneous assumptions, important 
local droplet concentration effects, and the necessity of 
incorporating film modeling. 

Current and future work is focused on: 1) improvement 
and validation of two fluid models (e.g., splashing), 2) larger 
time accurate models of shrouds with realistic penetrations, 
3) validation of multiphase analysis with forthcoming data 
from NASA Glenn’s Gearbox Windage Facility, 4) further 
evolution of film modeling, and, 5) incorporation of gear 
meshing dynamics and thermodynamics. 
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Summary 


As of the writing of this final report, our group’s research into gearbox windage aerodynamics 
continues. The physics understanding, and spoiler modification that it led to, have resulted in a US 
Patent. This patent has interested two major US rotorcraft manufacturers, Bell and Boeing, to invest in 
sponsored research of our group to explore its applicability to helical and spiral-bevel systems. Both of 
these efforts led to CFD demonstration that such modifications should improve the loss performance of 
rotorcraft transmissions. The exposure of our NASA sponsored research in this area to the wider 
aerospace community has led directly to sponsorship from Pratt and Whitney for exploration of the 
aerodynamics in their fan-drive gearbox. The multiphase modeling research supported in this NRA has 
similarities to some of the requirements of rotorcraft (and fixed wing) icing modeling (specifically, 
droplet and film physics), and this has benefitted our separate research in icing CFD modeling (to be 
published in this summer’s Asian Rotorcraft Conference). Lastly, our progress has led directly to our 
current 5 year VLRCOE project on rotorcraft loss-of-lubrication modeling, of which windage is an 
important component. 
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Executive Summary 

This document represents the evolving formal documentation of the NPHASE-PSU 
computer code. Version 3.15 is being delivered along with the software to NASA in 2013. 

Significant upgrades to the NPHASE-PSU have been made since the first delivery 
of draft documentation to DARPA and USNRC in 2006. These include a much lighter, 
faster and memory efficient face based front end, support for arbitrary polyhedra in front 
end, flow-solver and back-end, a generalized homogeneous multiphase capability, and 
several two-fluid modelling and algorithmic elements. Specific capability installed for the 
NASA Gearbox Windage Aerodynamics NRA are included in this version: • Hybrid 
Immersed Overset Boundary Method (HOIBM) [Noack et. al (2009)1 ■ Periodic boundary 
conditions for multiple frames of reference, ■ Fully generalized immersed boundary 
method, • Fully generalized conjugate heat transfer, • Droplet deposition, bouncing, 
splashing models, and, ■ Film transport and breakup. 

Overview of NPHASE-PSU 

NPHASE is a CFD code developed by Robert Kunz at the Penn State University 
Applied Research Laboratory (PSU-ARL) and Steve Antal at Rensselaer Polytechnic 
Institute (RPI). The code has been under development since 1998. Since NPHASE Version 
2.0 was established in 2000, two separate versions of the code have been developed 
independently by Kunz and Antal. This document applies to the version developed by 
Kunz at ARL Penn State, NPHASE-PSU which is distributed for free to three sponsoring 
US government agencies: NASA, DARPA and the USNRC. 

NPHASE-PSU is not a commercial CFD code, nor is it used for commercial 
consulting. NPHASE-PSU is open source. The mission of the software developer is to 
support government and industrial sponsors of programs related to PSU-ARL’s core 
research activities. The software is owned by the Pennsylvania State University and is 
distributed freely to research sponsors (including all source, tutorials, documentation, and 
front/back-end processing tools). 

NPHASE-PSU is written in standard ANSI-C, with a smattering of C++, and 
compiles under (at least) the open-source GNU C compiler, gcc/g++ and (in 2013) the intel 
LINUX compiler, ice. NPHASE-PSU refers to the CFD code itself, but employs several 
front-end and back-end processing tools for domain decomposition and reassembly, grid 
readers for standard COTS formats, pointer topology construction, and writers to standard 
postprocessing software file formats. These processing codes are written in FORTRAN 
77/90, ANSI-C, and C++. Accordingly, if the user wishes to modify these front/back-end 
programs they must also have access to a FORTRAN 90 compiler. NPHASE-PSU also 
requires several open source software libraries including MPI2, PETSC, METIS and 
SUGGAR/DirtLib each of which must be installed with NPHASE-PSU on the system. 
Although the code has in the past been installed on Windows and SGI systems, the present 
delivered version, V3.15, is verified to install and run only on desktops and clusters 
running LINUX. 
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NPHASE-PSU has the following characteristics, features, and capabilities: 

• Arbitrary number of fields and/or species, where different species are assumed to be 
in dynamic and thennodynamic equilibrium, and different fields are not (i.e. have 
different velocities and enthalpies). Mass fraction and volume fraction transport 
options are available for species/field transport. 

• Numerous interfacial mass, momentum, energy and turbulence exchange models 
associated with multiphase flow simulations. 

• 3D unstructured: Arbitrary polyhedral fonnulation with front-back ends supporting 
4 standard element types: tetrahedral, hexahedra, pyramids, prisms, as well as 
completely arbirtrary element types (n-faced polyhera). 

• Overset mesh capability, utilizing open source Suggar and DirtLib software. 

• Moving and deforming mesh capability (Geometric Conservation Law satisfying). 

• Fully matrix level parallelized using MPI2 and domain decomposition. 

• METIS used for domain decomposition embedded in front end. 

• PETSC and simple point linear equation solvers. 

• All-Mach number fonnulation: incompressible, weakly compressible, strongly 
compressible flows. Isothennal, Boussinesq and perfect gas single-phase 
compressible state relations are available. 

• Segregated pressure based algorithm and CPE algorithm for multiphase flow. 

• Face based finite volume scheme: 1 st through 3 ld order accurate convection 
discretization schemes, 2 nd order accurate viscous term discretization. 

• 1 st and 2 nd order, dual time based temporally accurate fonnulation. 

• Several low and high Reynolds number form 1 -equation, 2-equation, and v2f 
turbulence models. 

• Structural mechanics coupling to NASTRAN. 

• Radiation heat transfer coupling to RADTHERM. 

• Numerous “specialty” face and volume elements (conducting solid regions, porous 
regions, various quasi- ID conjugate heat transfer boundaries). 

• Full turbomachinery capability (rotating and non-rotating reference frames) 
including rotor-stator interaction and body force modeling. 

• “Light” face based file formats supported in front end. 

• ENSIGHT file format supported in back end. 

• Coded purely in ANSI-C, and C++ with some front and back end utilities coded in 
F77, F90. 

Partial development (features that are not fully implemented but are in source code 
in various stages of completion): 

• Non-isotropic mesh adaption 

• Full Reynolds Stress modeling 

• Conformation tensor transport 

• VOF for discrete interfaces 

• 6DOF dynamics 

• Fully coupled parallel algorithm 
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NPHASE-PSU has been applied to and validated against a broad range of complex 
single-phase and multiphase configurations including: 

• Gas-particle flows through a branching pipe junctions and human lung geometries 

• Bubble column reactors 

• Full-annulus rotor-stator pump and turbine stage analyses, including rotor-stator 
interactions 

• High Reynolds number submarine configurations at a range of angles of attack 

• Power plant cooling ponds 

• Microbubble drag reduction applications of relevance to DARPA 

• Geometrically complex UUV (MRUUV) and SEAL deliverivery vehicles (ASDS) 

• Several surface ship configurations (5415, Athena) 

• High speed maritime lifting pod 

• Micro-flows of biological cell systems 

• Numerous multiphase flows of relevance to the NRC (thermally driven counter- 
current reactor flows, 2-phase duct and pipe flows) 

• DES simulations of urban/atmospheric dispersion 

• Bubbly surface ship wakes 

• Thennal management of ta nk engine compartment 

• Thennal management of eco-friendly structures. 

• Droplet impingement in compressible aerodynamics flows for icing simulation 

• Gear aerodynamics 

Documentation of many of these cases appears in Kunz et al. (2001, 2003, 2007, 
2011) or can be obtained from the author. 
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NPHASE-PSU Theory Manual 


Governing Equations 

The single-pressure ensemble averaged continuity and momentum equations are 
cast in conservation law form as: 


d ak p k + d a P u j _ y /pik _ pki i 
dt ti ’ 


da k p k u k da k p k u k u k k dp 
— — -i — = -a 1- ■ 


at 


ax- 


aX- fiX; 


ot k pf 


auf au. 


kb 


■ + - 

ax- ax- 
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p k a k g, + M“ + X(D kl [ul- u k ]+ r lk u| - r kl u k ) 


k*l 


( 1 ) 

( 2 ) 


Superscripts k and 1 designate donor and receptor fields for mass transfer ( r kl ), and drag 
( D kl ) and non-drag ( M kl ) interfacial forces. In general each field, k, will have a different 
density, volume fraction, velocity and viscosity. For single phase flow, equations (l)-(2) 
reduce to: 


gp , 5 P u j _ (J 
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at a Xj axj dx- 

For homogeneous multiphase flow, it is assumed that the fields are in dynamic and 
thermodynamic equilibrium, and equations (l)-(2) reduce to: 


aa k p k + aa k p k u k 
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(5) 


( 6 ) 


j i j 

where the set of momentum equations is reduced to a single equation for the mixture. 
Superscript m represents mixture quantities, hi equations (l)-(6), a high Reynolds number 
fonn viscous term is assumed with dilatation and turbulence energy terms neglected 
(although these terms are available in NPFIASE-PSU). Energy and turbulence equations are 
considered below. 


Physical Modeling 
Generalize Field Transport 

The generalized n- field formulation in equations (l)-(2) can be applied to non- 
equilibrium multiphase flows in two ways. The more fundamental approach involves 
solving mass and momentum equations for each field that is present. For example, in the 
context of disperse bubbly flows, one could solve a single continuous liquid field and a 
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number of disperse fields, “binned” by size. In this approach each bubble field exchanges 
momentum with the continuous field through drag and non-drag interfacial forces which 
depend in magnitude on the local interfacial area density of that field, Ai nt =6a gas /Db (for 
spherical bubbles). This approach was used in our earlier work [ Kunz et ah, (2003, 2007)1 , 
where up to 1 1 bubble fields were solved. 

Interfacial Area Density Transport 

An alternative is to solve a single mass and momentum equation for each phase 
that is present and to accommodate the variation in dynamics due to phase interface 
evolution by modelling and solving for interfacial area transport. For example, in the 
context of disperse bubbly flows, a single gas field continuity and momentum equation 
would be solved, and an interfacial area density transport (IADT) equation would also be 
solved to determine a local mean characteristic diameter for the bubbles. This approach 
significantly reduces the model’s CPU requirements compared to solving an (N+l)-field 
system (N bubble fields). The numerical complexity associated with interfield transfer 
terms is also reduced considerably. 

Since mass transfer can be fully accommodated in the context of IADT (details 
presented below), the physical appropriateness of employing IADT rests on whether the 
interface dynamics can be sufficiently captured using a single local mean inter-phase 
interfacial area, with an assumed/modeled distribution of characteristic size/shape about 
that mean. This is demonstrated to be the case for an example calculation below. Currently 
in NPHASE-PSU, a generalized IADT formulation is available, with physical models in 
place to accommodate disperse bubbly flows related to Microbubble Drag Reduction 
(hereafter MBDR). In this context mass transfer corresponds to coalescence and breakup 
(between bubbles of different sizes). The IADT formulation in NPHASE-PSU is presented 
here, in that context, although any interface evolution (e.g., annular flow, droplet laden gas 
flows) can be modeled through addition of subroutines corresponding to those available for 
the disperse bubbly flow models currently available in V3.1. 

Following Hibiki et ah (2001), the IADT equation with source terms for breakup 
and coalescence can be written: 


where a; is the interfacial area density, u gj j are the gas phase velocity components, and Ob 
and O c are source terms for breakup and coalescence, respectively. The interfacial area 
density is defined as: 



dt dx : 


(V) 


6a g 


( 8 ) 


where a g is the volume fraction of the gas phase and D is the mean bubble diameter. The 
source terms are rates of change of interfacial area concentration, written as: 



( 9 ) 


8 
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where cj)b and (j) c are the rates of change of bubble number density (1/m s) due to breakup 
and coalescence, respectively. The factor \|/ depends on the bubble shape, here taken as 
spherical, so \|/=1/(36tc). The particular models used for breakup and coalescence for 
MBDR are presented below. 

Figure 1 illustrates that the dynamics of MBDR can be sufficiently captured using a 
single local mean gas-liquid interfacial area, with an assumed/modeled distribution of 
bubble size about that mean. Three MBDR cases are considered, corresponding to three gas 
injection rates at injector plates near the leading edge of a very high Reynolds number flat 
plate flow (see “HIP LATE” tutorial below). First, each case was run with three bubble 
fields using an approximation to the experimentally measured bubble size distribution. 

Then each case was run using a single gas field and interfacial area density as described 
above. For these comparisons no coalescence or breakup was incorporated so as to isolate 
the effect of the different interfacial dynamics modeling approaches. Details of the 
HIP LATE simulations are provided below, but Figure 1 serves to illustrate that 
incorporating interfacial area density has only a small impact on accuracy of drag reduction 
and bubble velocity predictions for MBDR. 



Q a /[Q a+ U e b0 o ] 


Figure 1. Comparison of 2-field and 4-field simulations for Uoo=18 m/s HIPLATE cases, (top) Drag 
reduction vs. x, (bottom) Normalized bubble velocity vs. normalized flow rate. 
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Interface Dynamics 
Overview 

The structure of NPHASE-PSU supports arbitrary forms for drag ( D kl ) and non- 
drag ( M kl ) interfacial dynamics models that appear in equation (2). The focus of two-fluid 
NPHASE-PSU research performed to date has been in the context of disperse bubbly flows 
(where the single continuous field is liquid) and disperse particle flows (where the single 
continuous field is gaseous). Accordingly the physical model set that currently resides 
within V3.1 are appropriate for these interface dynamics. 

A suite of bubble dynamics models have been developed, adapted from the open 
literature, and calibrated over the course of the development of NPHASE-PSU. These 
interfacial force models are summarized here. These models have been used in the context 
of a full-up n-bubble-field formulation where the bubble diameters that appear and hence 
the implied interfacial area between that field and the liquid are unique to and 
representative of that field. As indicated in the previous section, these models are also 
implemented in the context of a single gas field represented by a mean bubble diameter and 
attendant interfacial area density, which is transported and evolved with the flow. 

Drag 

In the context of particles, a conventional corrected Stokes drag law is incorporate 


D 


kl 


ip‘“ c o 


u!-u k 


6a gas _ 24 

It'D ~ 


D. 


Re. 


foK.) 


( 10 ) 


where the local particle Reynolds number is Re p = p gas 
drag-law correction used [ Loth (2000), for example] is: 
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rel 


D / p gas . The solid particle 


f D =l. + 0.1875Re p 
f D =1. + 0.1935 Re p 6305 


f D =1. + 0.015 Re p + 0.2283 Re p 
f n =0.44Re„/24. 


.427 


for Re p < 1 
for 1 < Re p < 285 
for 285 < Re p < 2000 
for 2000 < Re p <3.5xl0 5 


( 11 ) 


L D — p. 

In the context of bubbles, drag models have been implemented for spherical bubbles 
in seawater, clean fresh water and contaminated (tap) water. Again, a corrected Stokes drag 
law is employed: 


D kl = — p liq C D 
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u! -u k 


6a liq 24 f ^ 

a o a , =— ,C D =— f D (Re b ) 


D v 
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(12) 


where the local bubble Reynolds number is Re b = p 


- A4 


V 


rel 


D b /p m . 


For fresh water without impurities, the drag-law correction [ Loth (2000) , for 
example] is: 

f D =l. + 0.1875Re b forRe b <0.1 

f D =1. + 0.0565 Ref 5 for 0.1 < Re b < 500 H ' ' 

For contaminated (tap) water, the drag-law correction for solid spheres equation 
(1 1), is used. For seawater, a drag-law correction due to Detch (1991) is available. 


NASA/CR— 2014-2 1665 1 


10 

52 


In addition to water purity, locally high gas volume fraction and bubble deformation 
can influence the drag, so corrections to the spherical bubble, disperse flow models in 
equations (11) and (13) may be appropriate. For uniformly disperse flows, an increased 
drag coefficient is appropriate [Richardson-Zaki (1954), for example], and for flows where 
gas structures are streamlined (bubble columns, sheets) a reduced drag coefficient is 
appropriate. This latter effect likely is important in the near injector region of MBDR 
flows, where application of the standard disperse flow model gives rise to too much local 
drag, thereby inhibiting the penetration of the injected gas into the boundary layer. This 
observation became clear in the course of the FIIPLATE validation studies, where a 
significant defect in measured bubble velocity could not be obtained unless a “cluster” drag 
model was incorporated. Specifically, a model proposed by Johansen and Boysan (1988) 
has been adapted to an Eulerian framework: 


C n =C D0 | 1-1.54 


MIN (. 5157, a gas )] 23 


(14) 


where Cdo is the original drag coefficient in equations (1 1) or (13), a gas is the total gas 
volume fraction and the MIN function is provided to ensure that the corrected drag 
coefficient does not drop to below 1% of the uncorrected value. The importance of 
incorporating such a cluster drag fonn is demonstrated in Kunz et al. (2007) . 


For air-water droplet flows of aerodynamic interest, a droplet drag model is 
implemented Michaelides (2006) is employed: 


C d =0.4+ 


24 

Re, 


6 
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Virtual Mass 

Virtual mass is modeled following Lahey and Drew (2000) : 
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liq-gas _ gas 
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DV 
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liq 


(15) 


Lift 

The lift model employed in the NPHASE-PSU also follows Lahey and Drew 

( 2000 ) : 


Mi" = a gas p liq C, V rel x V x V liq 

LIFT L 1- - (16) 

Wall Lift 

An empirical turbulent near-wall bubble lift force has been implemented based on 
the formulation of Kawamura and Yoshiba (2004) . This force can be thought of as a 
repulsive force due to wall collisions. The form of the wall-lift force used is: 

Fwl = Cwl^ / 6)p' lq (k/D b )F damp 
F dam P = °- 5 [l - tanh(y waU / D b - 1 .5)] 

where Fdamp decays the force to zero away from the wall and the model constant used here, 
C WL =0.012/^1 + st k , is significantly smaller than that proposed by Kawamura and 
Yoshiba. The Stokes number is defined as St k =D B 2 p liq s/(18kp m ). 
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Turbulence Dispersion 

The homogeneous turbulence dispersion model is implemented in the framework of 
the Carrica, et al. (1999) gradient diffusion force model. The general expression for the 
dispersive force per unit volume (N/m ) may be written as: 


M 


k,TD _ liq V t r da 

1 “ P Sc D dx : 


(18) 


where C TD is the turbulent dispersion coefficient (units s' 1 ). For the Carrica, et al. (1999) 


•'TD 

model, C TD is defined as: 


C td - 


3 Cr 


8 R v 


Mrd 


(19) 


where uf el is the relative velocity between the continuous phase and disperse phase “k”. 


At the high gas volume fractions, dispersion is enhanced by collisions among 
bubbles. A new dispersion model has been developed, based on the collision frequency 
from the Prince-Blanch (1990) coalescence model. This dispersion mechanism is used in 
addition to one of the homogeneous turbulence dispersion models discussed above. Since 
DNS computations [Maxey, et al. (2005)1 show a significant effect of collision on 
dispersion for high gas volume fractions, a heuristic dispersion model based on the bubble 
collision rate has been formulated. The collision-induced dispersion model is implemented 
in the framework of the [ Carrica, et al. (1999)1 gradient diffusion force model, equation 
(18). 

We assume the dispersion model coefficient is an unknown function of the 
“dispersive collision rate”, which excludes bubbles that coalesce. To properly formulate the 
coefficient relationship, the collision rate must be normalized. For that purpose we choose a 
turbulent characteristic bubble response time ( x J BC ) defined as: 


x 


j _ 
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4 D: 


3 C 
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( 20 ) 


where k is the turbulent kinetic energy. Note that this is the bubble response time normally 
used to define the Stokes number, 


St 



( 21 ) 


Alternative characteristic times were evaluated with some success; however, the 
above relation is a reasonable choice with physical basis. 


The normalized dispersive collision rate (0™) for bubbles “i” and “j” with an 
equivalent volume V- is written as: 



\ 


0ij =0,jVij(l 
/ 

_ ^ij ^BC 
v 

(22) 

V = V: + V: 

1/2 


!J V 1 J 

/ 



NASA/CR— 2014-2 1665 1 


12 

54 


The turbulent dispersion coefficient (for a bubble “j”) is chosen to be 
proportional to the square root of the dispersive collision rate (nonnalized by a 
representative time scale, i J BC ). 


C 


coll 

TD 



(23) 


where C TD is a constant to be detennined. Note that the square root is chosen to obtain a 
consistent relation with the collisional pressure identified by Maxey, et al. (2005) . The 
Brown DNS calculations confirmed the functional dependence of collisional pressure on 
volume fraction. 


Further modifications to the dispersion model are required to treat other conditions, 
especially limiting cases with high gas volume fraction. A heuristic model as been 
implemented for the dispersion models and bubble lift model. 


The dispersion in NPHASE-PSU is modeled by summing the two contributions 
discussed above (equations (19) and (23)), i.e., 



In general this relation applies to each bubble field “k”. 


Deposition 


Deposition of droplets or particles can be modelled as inter-field exchange through 
mass transfer in the context of full-two fluid modelling or interfacial area density. However, 
in some applications, including particle deposition within the lung or droplet deposition onto 
an airfoil, it may be suitable to model the process as a mass sink. Specifically, once can 
assume that when the mass deposits on the wall it “leaves the domain” as in a transpiration 
condition, hi NPHASE-PSU this is handled using deposition models that simply extract 
mass and attendant momentum from the wall adjacent cells. 

For droplet deposition a simple inviscid impaction model is adopted wherein, all 
droplet momentum oriented into the wall is extracted along with its mass: 


^deposition = } P^Y^dA,- 

Awal/ 

where the droplet velocity is taken as the nodal value adjacent to the wall. This flux is 
consistently extracted from the droplet continuity and momentum equations. Other models 
are available in the code including diffusion based deposition for co-flowing systems, but 
this is not yet documented here. 
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Breakup and Coalescence - Multi-Bubble-Field Formulation 
Overview 

A general fonnulation for discrete bubble size distributions based on the approach 
of Kumar and Ramkrishna (1996 ) has been implemented in NPHASE-PSU. The approach 
allows one to rigorously conserve two functions of the bubble distribution function (or 
kernel) regardless of the discrete bubble sizes (bins) selected. There is a unique formulation 
for coalescence and breakup of bubbles. In both cases we have chosen to conserve volume 
moments of the bubble number density distribution function, n(v,t), i.e., 

oo 

M, = jVn(v,t)dv ( 25) 

0 

where v denotes the bubble volume and t is time. Of course, the distribution function is a 
function of spatial location as well. We have chosen the zero-th (p=0) and first (p=l) 
moments at present, though the coding permits arbitrary moments to be conserved. The 
rational for this choice is the conservation of the number of bubbles and the volume of 
bubbles during the coalescence and breakup process. Though the interfacial area of the 
bubbles is an important quantity in two-phase bubbly flows, it is not conserved in general. 
The correct interfacial area will be preserved by conserving the number of bubbles and the 
volume of the bubbles. 

It should be noted that other investigators [e.g., Carrica, et al (1999)1 have used 
fonnulations based on bubble mass, since in cases with significant gas compressibility the 
bubble mass is conserved while the volume changes (in the absence of either coalescence 
or breakup). The present implementation of the Kumar-Ram kr ishna scheme in NPHASE- 
PSU easily permits the use of bubble mass as the bubble size metric rather than bubble 
volume, if necessary. 

Prince-Blanch Coalescence Model 

For coalescence, the rate kernel employed is due to Prince and Blanch (1990) and 
Williams and Loyalka (1991) . The latter text offers a fairly complete description of the 
physics of coalescence and various mathematical approaches for modeling the various 
coalescence mechanisms. Three primary mechanisms may be included in the coalescence 
kernel - (1) turbulent diffusion, (2) “laminar” shear, and (3) buoyancy. The so-called 
“laminar” shear contribution is modeled as a function of the local velocity gradient, and is 
relevant only for laminar flow and therefore, is not considered. The fonnulation of Prince 
and Blanch models the effect of turbulent diffusion due to “small” eddies, while the 
fonnulation of Williams and Loyalka also purports to model the effect of small eddies, 
though with an approach that relies on a bubble scale that is small compared to the scale of 
the turbulence. Hence the Williams and Loyalka formulation may not apply to the bubble 
sizes expected to be present in microbubble drag reduction applications. 

The turbulent diffusion contribution is due to a statistical average of the fluid 
velocity fluctuations. However, a general turbulent flow also has a mean velocity gradient 
which has an effect on collisions. Williams and Loyalka discuss the impact of a laminar 
flow velocity gradient on coalescence. For the present application their formulation was 
adapted to treat the mean velocity gradient effect in turbulent flow. 

The coalescence model considering turbulent diffusion due to small-scale 
turbulence and mean-shear is operational in NPHASE-PSU. The effect of buoyancy on 
coalescence has been neglected. 
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The turbulent collision rate is a dominant factor in bubble coalescence according to 
both Prince and Blanch (1990), and Williams and Loyalka (1991) . For small eddies the 
turbulence is assumed to be isotropic (at least on the scale of the bubble diameter) and the 
bubble size is assumed to lie in the inertial subrange. The same assumption is made in the 
breakup model fonnulation discussed below. Following Prince and Blanch, the collision 
frequency [ <9 ;/ ,1 /( nr s)] between bubbles i and j due to turbulent motion may be written: 


/ — — \ 1/2 

e^n^S^uf+uj) (26) 

where n, and n j are the number densities (nT ) of bubbles with diameters D, and Dj, 
respectively. Also, u] is the root mean square of the fluctuating velocity of bubble i and 
Sy is the collision cross-sectional area defined by Prince and Blanch: 


s»^( D , + D ,) 2 

(27) 

The required fluctuating velocity in the inertial subrange for isotropic turbulence is 
given by Prince and Blanch: 

u 1 = V2(sD i ) 1/3 Ui =V2(eD i ) 1/3 

(28) 

where the relevant turbulence length scale is assumed to be the bubble diameter. T 

re 


leading constant in equation (28) is not universally agreed upon in the literature and the 
turbulence length scale also appears in several different forms, although always as a 
function of bubble diameter. 


Combining the above expressions yields the desired relation for the collision 
frequency, 


T 

Vffii 

16 

-s 1/3 (d. +D.) 2 (D 2/3 +D 2/3 ) 1/2 

(29) 


The probability that a collision results in coalescence is required to complete the 


rate kernel formulation. Again, following Prince and Blanch, this probability is tenned the 
collision efficiency and is a function of the contact time between bubbles and the time 
required for bubbles to coalesce. For a pair of bubbles, this efficiency ( X- ) is written as 
[following Coulaloglou and Tavlarides (1977)1 : 


k y =exp(-t ij /T y ) 

(30) 

where A is the time required for bubbles of diameters D, and Dj to coalesce and f (/ . is the 
contact time for the two bubbles. From other literature, Prince-Blanch presented the 
following expression for the coalescence time ( t Lj ): 

t. = 1 J ' In 0 

1J 16a l h f J 

V J 

(31) 
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where h 0 is an initial film thickness between two bubbles as they just come into contact 
and h , is a final critical film thickness where rupture occurs and the bubbles coalesce. The 

quantity D tj is an equivalent diameter for bubbles of unequal size and is given by: 


D u = 

f 2 2 ^ 

iL + cJ 

-1 

(32) 

For air-water systems, the film thickness values quoted by Prince-Blanch (j 
other sources) are 

'rom 

h 0 =l(T 4 m,h f =l(T 8 m 

(33) 


Finally, an estimate of the contact time ( Ty ) for bubbles in turbulent flow was 

made by Levich (1962) from dimensional analysis. A modification due to the relative 
velocity between the bubbles is noted by Carrica, et al. (1999) , resulting in the following 
expression: 



where D ch is a characteristic length related to the bubble sizes and u re ijj is the mean 


relative velocity between the colliding bubbles. The characteristic length ( D ch ) in equation 
(34) may be taken as an adjustable parameter in this model, hi the absence of better 
information, D ch will be taken as the average of the inverse of the bubble diameters, 

2 D ch = ( IT 1 + Dj 1 ) ' , as suggested by Carrica, et al. (1999) , and Prince and Blanch (1990) . 


Furthermore, all quantities in the model are assumed to be statistical averages for a 
turbulent flow, thus further uncertainties in the model may result. There appears to be very 
little data or analysis in the literature addressing these complex issues. 

Lehr-Mewes Breakup Model 

For breakup, one rate kernel investigated is due to Lehr and Mewes (2001) . This 
kernel has some important properties that allow the formation of a small bubble and a large 
bubble when a large bubble breaks up. The breakup mechanism considered is due to small- 
scale turbulent eddies. The Kumar-Ramkrishna (1996) fonnulation for breakup requires 
evaluation of volume integrals of the rate kernel, and the fonn of the kernel has some 
characteristics that can lead to numerical problems if not addressed carefully. 

The Lehr-Mewes rate kernel for the (binary) breakup of a “mother” bubble with 
non-dimensional volume x k into daughter bubbles with non-dimensional volumes v and 

(x k - v) is given by: 


x 1/3 

r i(v,x k ) = C J /3 

F . (v) 1 

mm V v / „ 7/9 

X k 

for v < — 
2 

(35) 


NASA/CR— 2014-2 1665 1 


16 

58 


The non-dimensional daughter bubble volume v is a defined as: 


„ V 

V s, 

(36) 

and v st is related to the maximum stable bubble, v stable , size by: 

tt rr 9 / 5 

v = 2 3/5 = 2 3/5 V 

stable 6 p 9/5 e 6/5 st 

(37) 

where a is the surface tension (N/m) between the gas and liquid phases, p , is the liquid 

density (kg/m 3 ) and £ is the turbulence energy dissipation rate (m 2 /s 3 ). The function F min 
is given by: 

Wv) = v 7 ' 6 for v < 1 

(38) 

Fmin(v )-„ 7/9 for V > 1 



Note that the rate kernel is symmetric about v = x k / 2 which allows its evaluation 
for v > x k /2 using equation (35). Since the rate kernel must be non-negative, equation 
(35) must be restricted. This leads to a minimum value for a minimum daughter bubble size 
given by: 


v - x~ 2/3 

V min X k 


(39) 


where the rate kernel becomes zero. Also, the rate kernel has a slope discontinuity at 


v = v_ 


The fonn of the function F m[n also gives rise to a slope discontinuity in the rate 
kernel at v = 1 , and possibly at v = x k /2. Furthermore, the kernel has very large gradients 

for large non-dimensional bubble sizes. This is shown in Figure 2, where the normalized 
daughter size distribution for the L-M rate kernel is shown for several values of a volume 
ratio parameter, VR, defined by: 


VR = ^_ 


Atable ^stable 


(40) 


where v slable is the maximum stable bubble volume. As a result of equation (37), the non- 


dimensional bubble volume is a function of the local flow properties, thus it will vary 
throughout the flow field. 


In general the required moment integrals of the rate kernel required for the K-R 
fonnulation cannot be evaluated analytically. The zero-th moment is an exception. 


All moments can be evaluated numerically; however for large values of VR, 
accuracy has been shown to be poor unless caution is exercised in selecting the integration 
step size. This is due to the large gradients shown in Figure 3. An approximate analytical 
evaluation of the kernel integrals for large VR was explored, but proved to be impractical 
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and did not reduce CPU time. Thus an adaptive procedure was implemented for selecting 
the integration step size based on a prescribed accuracy in resolving the L-M rate kernel. 
For a very wide range of mother bubble sizes, this approach requires only a moderate 
number of integration steps (< 1000) to detennine the necessary moments very accurately 
(<0.01% error). 

Martinez-Bazan Breakup Model 

Another rate kernel investigated is due to Martinez-Bazan, et ah ( 1999 a, b) . This 
kernel is much different than the Lehr-Mewes kernel in that the formation of a small bubble 
and a large bubble from a bubble breakup has very low probability. The breakup 
mechanism considered is due to turbulent eddies and a phenomenological model for the 
breakup kernel (frequency) was developed using experimental data from a high-Reynolds 
number water jet flow with bubble injection. The experiments were conducted very 
carefully to insure that the turbulence in the jet was locally homogeneous, isotropic and in 
near-equilibrium. The model assumes that the initial bubble size, D 0 is in the inertial 

subrange, i.e., rj « D 0 « L x , where 77 is the Kolmogorov microscale and L x is the 
integral scale of the turbulence. 



( 

X 


n = 

l) 


( 41 ) 

l £ J 
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x 


rc E iifri_=°) 
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( 42 ) 
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v/x k 

Figure 2. Normalized daughter size distribution for Lehr-Mewes rate kernel. 



v/x k 


Figure 3. Normalized daughter size distribution for Lehr-Mewes rate kernel near V/Xk=0, for VR=140. 
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where u' is the fluctuating component of axial velocity, k x is the turbulence wave number 
in the axial direction and the tensor E is the turbulence energy-spectrum function [ Hinze 
(1975)1 . 

It should be noted that the experimental technique of Martinez -Bazan, et al. (1999a, 
b] had a minimum measurable bubble size of 83 pm, which Martinez-Bazan states did not 
affect their breakup frequency results. 


A critical bubble diameter D c exists and if D <D C the bubbles will never breakup: 


( 

D c = 

V 

12a A 

Pp ) 

3/5 

s- 2/5 

(43) 

A minimum diameter exists below which there is insufficient turbulence induced 
stress to result in bubble breakup. 

D mm = 

f 12a ") 
IppD; 

3/2 

8 1 

(44) 

The breakup frequency (1/s) is given by: 

. Jp(sD) 2/3 -12o/(pD) 

g(s.D) = V V ’ D A ’ 

(45) 


where f> = 8.2 [Batchelor (1956)1 and K g = 0.25 was detennined experimentally by 
Martinez-Bazan, et al. (1999a) . 


The Martinez-Bazan, et al. (1999a) breakup model was developed for conditions 
more representative of MBDR flows than the Lehr-Mewes model, thus the fonner has been 
utilized in the present work. 

Kumar-Ramkrishna Partitioning-Breakup 

The Kumar-Ramkrishna particle bin size representation is shown schematically in 
Figure 4. Here x, is the representative bubble bin volume (e.g. average or mid-point 
volume) due to breakup of mother with volume x k 
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Figure 4. Kumar-Ramkrishna particle bin size representation. 


Conservation leads to the following relation with the breakup rate is written as: 


R 


Bb 


Z n i,k r B ( x k) N k(t) 


(46) 


NASA/CR— 2014-2 1665 1 


20 

62 


where N k is the total number of particles in bin “k” and V H (x k ) is the breakup frequency 
(kernel) for mother particle x k and n i k is the contribution to the population of the “z'-th” 
bin size ( x i ) due to breakup of particle x k . 



Bf k = jVp(v,x k )dv 

x i 

(48) 

The death rate due to breakup of a particle of size x k is: 

R Db =r ( X k) N k (t) 

(49) 

Kumar-Ramkrishna Partitioning-Coalescence 

The coalescence formulation is simpler than that for breakup. The birth rate of 
particles due to the coalescence of particles in bins j and k is given by: 

j>k 

R d = £ (i-o.58,. k )n J . l r' k N J N k 

j.k 

X M <(Xj+X k )<X i+ i 

(50) 

where the distribution function due to the coalescence ( r] jk ) is given by: 

V^ x v _v v x p 

u - i+1 i+1 X < V < X 

^ Y H Y V Y v M » X , SVSX i + l 

A i A i+l A i A i+l 

vX-rv’x^ ^ ^ 

T l J .k= „ v V X i-1 - V - X; 

x r x i-.- x i x w 
where v = Xj + x k 

(51) 

where v=Xj+x k , and T l,k is the coalescence rate (kernel) due to the coalescence of particles 
in bins j andk. 

The death rate of particles due to the coalescence of particles in bins j and k is given 
by: 

R cd = N iZr il N k 

k=l 

(52) 


Breakup and Coalescence - Interfacial Area Density Transport Formulation 

An approximate fonnulation including bubble breakup and coalescence within the 
interfacial area framework was proposed by Lehr and Mewes (2001) . Lehr and Mewes 
solved the population balance equation “to describe the evolution of bubble sizes in two- 
phase flow.” To reduce the numerical complexity due to a large number of equations and 
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strong coupling, they formulated an equation for average bubble volume (equivalent to the 
interfacial area transport equation) using an approximate analytical approach. A summary 
of the Lehr-Mewes approach follows. Source terms in the population balance equation 
involve breakup and coalescence kernel functions that are a function of the bubble volume, 
v. By assuming that an arithmetically averaged bubble volume ( v) may be used in the 
kernel functions, a simplified solution for the bubble number-density distribution function, 
f(v) , results: 


a ( y\ 

f(v) = = |exp 

V- ^ v) 

(53) 


00 (X 

n B = Jf(v')dv' = — 
0 v 

(54) 


Lehr and Mewes obtained a transport equation for average bubble volume with 


simplified source terms due to breakup and coalescence (equivalent to the source terms O b 
and ®c and in equation (7). The bubble number-density PDF implies a bubble size 
distribution consistent with the above noted assumptions. We use this PDF to evaluate 
bubble number densities for discrete “bins.” The bins are defined as shown in Figure 4. 


Here x; is the representative bubble bin volume (e.g. average or mid-point volume) 
of bin “i” and Vm and V; are the lower and upper bin volumes of bin “i”, respectively. The 
number density PDF of bubbles in bin “i” is then 


N =— (e~ 


Vj_] / V 


-v,/v) 


(55) 


This result approaches the number density PDF for a sufficiently large number of 
bins and a sufficiently large maximum bin volume. Also the first bin is assumed to contain 
all bubbles from zero bin volume to the uppennost volume of this bin (i.e. Vi = 0). Further 
to prevent errors due to an insufficiently “large” maximum volume, the distribution must 


normalized such that V m = i . 

Zj B(i) 1 
all bins 


As in the N-bin formulation, we use the Prince and Blanch (1990) rate kernel for 
coalescence and the Martinez, et al ( 1999a, b) rate kernel and daughter size distribution for 
breakup. A complete description of these models is included above; only the essentials are 
summarized here under the assumption that the rates may be evaluated using the mean 
bubble diameter. 


4>c = n 


UJ 


1 / 3 — 7/3 

j 1/3 D 


ex 


pHb/^b) 


(56) 


where n B is the bubble number density, e is the turbulence energy dissipation rate, t B is 


the time required for two bubbles of diameter D to coalesce and T B is the contact time for 
the two bubbles. In the interfacial area density formulation, the bubble number density is 
given by 
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n B- -2 

7lD 

(57) 

As in the prior section the time required for two bubbles to coalesce is given by: 

t B = 

( / — \3 ,. ^ 

(0.5DJ p liq 
16a 

1 ) 

1/2 

l h f J 

(58) 


where h a is an initial film thickness between two bubbles as they just come into contact 
and h f is a final critical film thickness where rupture occurs and the bubbles coalesce. For 

air-water systems, the film thic kn ess values quoted by Prince and Blanch (from other 
sources) are 


The contact time for bubbles in turbulent flow [ Levich (1962)1 with a modification 
due to the relative velocity between the bubbles [ Carrica, et al. (1999)1 , is given by: 


D 


X B 


ch 


U B ,„,+2(0.5D d ,e) l ' J 


(59) 


where, as before, D ch is a characteristic length related to the bubble sizes and u Brel is the 
mean relative velocity between the colliding bubbles. The characteristic length ( D ch ) is 
taken as D ch = D . 


Enthalpy Transport 


For compressible flows and flows with heat transfer, it is necessary to solve for an 
energy equation. NPHASE-PSU incorporates an enthalpy transport equation for each field: 


5 (a k p k h k )+ 5 (a k p k u k h k ) = a kDLp + 5 
5t v ’ 5x j v J ’ Dt dx j 

" J k p k W~ 

or [i k +7h_ 

l aJs Xj 

kSuf 
+ Ty dx } 

(60) 


An available alternative that is better suited for higher Mach number compressible 
flows (transonic, supersonic) in the context of the segregated solution strategy employed in 
NPFIASE-PSU is stagnation enthalpy transport: 


— (a k p k H k ) +— (a k P k u k H k ) = a k — + — 

a k 


1 

+— (x k u k ) 

(61) 
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Turbulence Model 

NPHASE-PSU has a number of low and high Reynolds number fonn 1 -equation 
turbulence models (Spalart-Allmaras), 2-equation turbulence models (k-s, q-co, k-co, k-R) 
and a low Reynolds number 4-equation v2f model [ Durbin, (1991)1 . In the context of 
multifield flows, separate turbulence transport scalars are solved for each field. For 
example, the high Reynolds number k-s model is written: 
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(62) 


In equation (62), all field indicator superscripts are eliminated if only the liquid 
field is solved. Sk and S E are available source/sink terms to: extract turbulence energy 
associated with breakup HVleng and Uhlman (1998), Kunz et al. (2003)1 , and modified 
production due to interface dynamics and mass transfer mechanisms proposed by various 
authors rFerrante and Elghobashi (2004, 2005) , Tryggvason and Lu (2005)1 . 

Numerics/Code 

For single phase flow, the present algorithm follows established segregated pressure 
based methodology. A colocated variable arrangement is used and a lagged coefficient 
linearization is applied [Clift and Forsyth ( 1994), for examplel . One of several diagonal 
dominance preserving, finite volume spatial discretization schemes is selected for the 
momentum and turbulence transport equations. Continuity is introduced through a pressure 
correction equation, based on the SIMPLE-C algorithm fVan Doonnal and Raithby 
(1984)1 . In constructing cell face fluxes, a momentum interpolation scheme I Rhic and 
Chow (1983)1 is employed which introduces damping in the continuity equation. At each 
iteration, the discrete momentum equations are solved approximately, followed by a more 
exact solution of the pressure correction equation. Turbulence scalar and volume fraction 
equations are then solved in succession. As discussed above, several important numerical 
issues arise in two-fluid CFD, foremost among these, that sufficient implicit coupling 
between the constituents be established. In the present work this is accomplished using the 
Coupled Phasic Exchange (CPE) algorithm Kunz et al. (1998) . In NPHASE-PSU, CPE 
has been extended to a fully unstructured, parallel, time accurate scheme employing higher- 
order discretization practices. Details of the data structure, discretization, and CPE 
elements of the scheme are summarized in this section. 

Data Structure 

The hierarchal data structure employed is illustrated in Figure 5. The cell-centered 
finite volume flow solver accepts arbitrary polyhedral elements. The data structure is face 
based, that is, subsequent to the assembly of geometric parameters in the front end, all 
inter-element connectivity is retained in face pointers to the two adjacent cells. The 
fundamental data structure member is the “fedge” (face edge) which points to its two 
vertices and faces. Each face points to its bounding elements. This data structure provides a 
convenient framework for assembly of all required geometric parameters. 

fedges and faces are identified as either internal or boundary. A “boundary_patch” 
structure in the C flow solver includes as members a number of attributes for boundary 
faces including areas and other geometric information, scalar values at the face center, 
fluxes, and inter-partition boundary data storage and transfer buffers. 
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Discretization 


The governing equations are discretized using a cell centered finite volume method 
applied to arbitrary polyhedral cell types. Inviscid and viscous fluxes are accumulated by 
sweeping through internal and boundary faces. For inviscid flux evaluation: 


|a k p k ()) k V k *dA = ^Cf(j)f 


(63) 


vertex 




Figure 5. Heirarchal data structure in NPHASE-PSU 

K k • k 

where Cj) is face mass flux for field k, and c|)| is the value of general transport scalar (j) 

evaluated at face, f. The summation is taken over all faces bounding the element. C ^ is 

evaluated based on field variables available prior to the solution of the transport equation 

for (j) k (lagged coefficient linearization). Second order accuracy is obtained by evaluating 

1 / 

Cj: using a central plus 4th difference pressure artificial dissipation term due to Rhie and 
Chow (1983) : C 



u 

and by evaluating cf> ^ from Lien (2000) : 
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=4>t + (V4» k -dr) u 

(65) 

In equation (64), the overbar denotes a geometrically weighted mean at the face, 
i.e., referring to Figure 6. : 

Vp = (l-s)(V Pl ) + s(Vp 2 ) 
s = 8s 1 /(8s 1 + 8s ? ) 

(66) 



Figure 6. Geometry nomenclature for cell face evaluations. 

and A designates a difference across the face (i.e., Ap=p 2 -pi). In equation (65) subscript U 
designates the quantity associated with the element upwind of face, f (which can vary with 
field), and dr is the vector from the upwind cell center to the face center. In Figure 7, 
results of a two-dimensional inviscid parallel stream test case are presented using a square 
mesh on a square domain aligned 45° skew to the flow direction and also using a triangular 
mesh. Inflow axial velocities are specified as =2 along the upper half inlet and =1 on the 
lower half. On both meshes, the significant interface smearing associated with first order 
upwinding is significantly reduced using the second order expression in equation (65). As 

detailed in Kunz et ah (1998) , dissipation parameters, B k and F k in equation (64) are 
scaled in a fashion that accommodates interfacial drag, mass transfer and dispersion forces: 



where P is the NPx NP point coefficient matrix for the momentum equations defined 
below, which incorporates drag and mass transfer. F k is consistent with a widely used class 
of dispersive interfacial forces, M kl = KVa 1 re.g., Lopez DeBertodano (1998)1 . 
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Figure 7. Comparison of first and second order convection discretization for an inviscid “mixing” 
layer. Flow is left to right. Axial velocity contours, white: V = 2 , black: V = 1 . 

The evaluation of B k and F k in equation (67), requires the inversion of a rank NP 
matrix P at each grid point, at each iteration. This potentially CPU intensive procedure is 
circumvented by applying a simple Jacobi fixed point iterative procedure to approximately 
invert P. This procedure is rapidly convergent (2 sweeps are employed) since the P 
matrices are very well conditioned as discussed below. 


Neglecting cross-diffusion and dilatation, the viscous flux in the momentum 
equations can be written for an element face as: 


( 68 ) 

Referring to Figure 6, the gradient of a scalar, cj), on the face can be written as: 
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The terms labelled A represent components of the gradient that are orthogonal to 
S 12 . These terms are generally small (for hexahedral or prismatic elements extruded from 
geometric surfaces, neglecting them is nearly equivalent to the thin-layer assumption). 
Their discrete form is treated explicitly in the solution of the momentum equations (tenn S k 
in equation (76) below). The terms labelled B represent components of the gradient that are 
parallel to s 12 . These are discretized as: 


J 

f 

(apVV'dA) = (ap) f 

( , 7 

vvA, 

l d *lJ 

( £ s *dA) - (ap) f ^ 

(^(ds'dA) 

ds 2 

(70) 

and are treated implicitly (terms A p and A k b in equation (76) below). 

Gradients that appear in the flux calculations, and elsewhere, 
Gauss’ Faw: 

are computed using 

V f 

(71) 


with internal face values of ()) f computed from equation (66), and the summation take over 
all faces bounding an element. Equation (71) is computed by sweeping all internal and 
boundary faces, accumulating adjacent element contributions to V and A f (j) k from the face. 


Interfacial Force Evaluation 


In order to discuss interfacial force discretization issues, we consider three classes 
of these terms. First, when cast as in equation (10), drag can be viewed as a scalar sink 
term. That is to say, drag term, D kl , which appears in the momentum equations as: 


ZD“(u|-uf) 


(72) 


is generally evaluated for each element, and by virtue of its relative velocity factor, 
incorporated implicitly in the NPx NP block diagonal, P, of the momentum equation 
coefficient matrix, as seen in equation (80) below. As discussed above, the appearance of 
drag in this term is accommodated consistently in the Rhie-Chow scale factors in equation 
(67). It is noted that numerically, mass transfer plays a very similar role to drag (though 
D kl = D lk and in general T kl ^ T lk ), and accordingly its treatment is consistent with that 
discussed for drag. 


The second class of interfacial force terms are those that are linear in the gradient of 
volume fraction. These are generally dispersive in nature. These terms are evaluated 
straightforwardly using model equations such as (18), however, as is demonstrated in Kunz 

et al. (1998) , including the Rhie-Chow-like term, p* a* F I Va*A f - Aa 


A 1 


in 


equation (67) is critical for obtaining convergent oscillation free solutions when such forces 
are present. 


The third class of forces are simply those that do not conform to drag-like or 
dispersive-like forms. An example is lift, a particular form of which is taken here from 
Fahey and Drew (2000) : 
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Mlift = a d p c C L y rel x V x v° 


(73) 


where superscripts c and d refer to continuous and disperse fields respectively. 


A straightforward discretization of equation (73) in an element centered (or variable 
colocated) scheme such as presented here, would involve evaluating gradients using 
equation (71) and multiplying by appropriate velocity, volume fraction and density factors 
using element values. In Kunz and Venkateswaran (2000) it was demonstrated that such an 
approach can also lead to solution oscillations and attendant convergence degradation. 
There it was observed that staggered grid methods (i.e., those where the momentum 
equations are evaluated at locations staggered to the element centers) do not exhibit this 
behavior for this class of force. Accordingly, a staggered force discretization was proposed 
wherein the force in equation (73) is evaluated at each cell-face. Face values are then 
averaged to obtain element values. This force distribution renders staggered and colocated 
forms identical for linear forces and thereby removes solution oscillations. In the present 
unstructured framework this “distribution” of force across several nodes can be written: 


I M ,V, 

M = — 

V 


(74) 


where Mf represents the force averaged to the face per equation (73), V is the element 
volume and V f is the volume formed by face f and the segments connecting the face ver- 


tices to the volume centroid. It was observed in Kunz and Venkateswaran (2000) that this 
approach is equivalent to the addition of a second difference artificial dissipation to the 
standard colocated discretization, i.e., in the present unstructured context: 


M = M + V*KVM 


(75) 


where scaling factor, K, has dimensions of length". 


Boundary Conditions 


A palette of boundary conditions are available in the code including walls, 
symmetry boundaries, inlets (transport scalars specified, pressure extrapolated from the 
domain interior), pressure boundaries (transport scalars extrapolated from the domain 
interior, pressure specified), and cyclic boundaries (for turbomachinery analysis). All 
boundary conditions are treated implicitly in the fonnation of influence coefficients for the 
transport scalars. For scintered metal plate injection, porous wall boundary conditions are 
used, where an area penneability, X, is specified. Shear force on porous boundary faces is 
apportioned as F= x w A|( 1 - X), where Af is the face area and A fk is the area available for 
injection flux. 


Implicit Solution Procedure 


Invoking a dual-time formulation, the discretized governing equations for transport 
scalar, cf) k , can be written in A-form as: 
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In equation (76), second order backward differencing has been used for the physical 
time derivative (At) and Euler implicit differencing is employed for the pseudo-time 
derivative (Ax). A standard under-relaxation procedure is employed where an appropriate 
underrelaxation factor, co is selected (0.3 < co < 0.7) and the pseudo-timestep is evaluated 
from: 



f > 



Ax- “ 

p k a k V 


(77) 

1 — CO 

A k + » kl 




V k*l 7 



It has been observed in Venkateswaran et al. (1997) that such a specification is 
equivalent to a local timestepping procedure that accommodates CFL and VonNeuman 
stability. For physical transients, pseudo-timesteps correspond to sub-iterations of the 
SIMPLE-C algorithm. 

Phase Coupled Scalar Linear Solution Strategy 


Equation (76) represents a coupled system of NP equations for the NP unknowns 

i k. 

<P : 


A(j) = S 

(78) 

where coefficient matrix A has the fonn: 
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and upper and lower block triangular matrices U and L containing neighbor cell influence 


coefficients. 


For the diagonal dominance preserving discretizations employed, conventional 
iterative schemes will have diagonally dominant iteration matrices with spectral radii less 
than or equal to the underrelaxation factor, co [Kunz et ah (1998)1 , a direct consequence of 
the well conditioned nature of the main diagonal block matrix P. Accordingly, we consis- 
tently employ a simple point Jacobi scheme for solving equation (76) for all scalars (u,, a, 
k, s), as this scheme is guaranteed to provide adequate convergence within several sweeps. 
For the momentum equations, all three velocity components are solved for all fields 
simultaneously using point Jacobi iteration. 

As discussed above, the well conditioned nature of P renders determination of 
dissipation parameters B k and F k , in equation (67), (which scale with P' 1 ) amenable to a 
simple point Jacobi iteration as well. 

Continuity Equation Linear Solution Strategy 

In the present work a mixture volume conservation equation is derived by summing 
individual field volume fraction equations, each normalized by field density. A SIMPLE-C 
[Van Doonnal and Raithby, (1984)1 based pressure-velocity corrector relation (which 
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accommodates the same level of interfield coupling as the artificial dissipation operators 
discussed above, Kunz et al. 119981) is applied to develop an elliptic pressure correction 
equation. Transport equations for the field volume fraction equations are then solved. In 
this fairly standard method, under-relaxation is not employed for the pressure corrector 
equation in order to achieve a measure of mixture volume conservation at each pseudo- 
timestep. As as result, the discrete pressure corrector equation system is symmetric positive 
semi-definite (A p =^A nb ) and thereby its linear solution is a challenging and important 

nb 

factor in the nonlinear convergence rate of the overall scheme. In NPHASE-PSU, the 
PETSC suite of solvers are employed for the solution of this system. Depending on the 
degree to which mass conservation needs to be satisfied at a given non-linear iteration, a 
GMRES solver or a more CPU intensive Algebraic Multigrid procedure are invoked from 
the PETSC library of solvers. Details of these solvers are provided in (PETSC 120061) . 

Parallelization 

The code is parallelized based on domain decomposition using MPI. Partitioning is 
carried out in the pre-processor, fump, as described in the user’s Manual below. Inter- 
partition boundaries are input to the flow code from fump as any other boundary condition 
with a single additional boundary patch attribute being the neighbor partition processor 
number, fump writes inter-partition face pointers to the NPHASE-PSU input files 
(unphase. gridxxx) in the same order that these faces are encountered in fump. Accordingly 
no reordering is required when loading and unloading 1-D structures associated with 
message passing. Data is passed after each scalar is computed in the segregated procedure. 
For the point iterative solvers used for the scalar equations, Acj) is passed at every sweep of 
the linear solver, so that there is no degradation in convergence due to domain 
decomposition. For the PETSC solvers used for the pressure corrector equation, the code is 
parallelized at the matrix level. Accordingly, a global matrix is assembled each non-linear 
iteration and global mass conservation is strictly enforced each timestep as the pressure 
solver is converged. 

Further details on the physical models, numerics and code are available in Kunz et 
al. (1999. 2000, 2001, 2007. 2011) . 
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NPHASE-PSU User’s Manual 


Preprocessing 

Overview 

Figure 8 illustrates the front end for NPHASE-PSU. The code accepts a “light” 
faced-based grid specification. This achieves several goals. Firstly it enables the 
specification of arbitrary polyhedral meshes rather than being restricted to the four 
simplicial element types (tetrahedral, hexahedra, prisms and pyramids). Secondly, most 
commercial grid generators produce the face-based COBALT files now native to 
NPHASE-PSU (Gridgen, Pointwise, ICEM, HARPOON) or closely related face based 
fonnats (e.g., GAMBIT). 
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Input Files 

Referring to Figure 8, there are three sets of input files to NPF1ASE-PSU. The first 
is the author file, nphase.dat, which is a simple key-word based ascii file that specifies all 
of the real and integer data defining execution control, boundary condition flags and 
attributes, fluid properties and initial conditions. A simple example nphase.dat file is 
included in Figure 9, and the “Control Commands” section below is devoted to a 
description of all of the key words available. The keywords (such as “number of fields”) 
are case sensitive but blanks are ignored. Integer and real attributes can be written in free 
format. Lines are commented out by placing # in column 1 . 


#case title: 

#simple nphase.dat file 

iterations to perform 100 

number of fields 1 

time accurate simulation 
temporal discretization momentum 1 
physical timestep in seconds .1 
number of physical timesteps 1 
transient file write frequency 2 

#initialize run with restart file 

produce ensight output 
restart file write frequency 100 

dont perform wall match logic 

inlet patch 1 0 

1 . 0 . 0 . 1 . 1 . 0 . 10 . 10 . 0.00 

pressure patch 1 0 

0 . 1 . 1 . 0.1 0.1 0 . 0 . 

turbulent flow high reynolds number k epsilon 

constant fluid molecular viscosity 1 . Oe-2 
constant fluid density 1.0 

function entry/exit echo off 

solversweepsf oru 3 
solversweepsf orv 3 
solversweepsf orw 3 
solversweepsf or k 3 
solversweepsf ore 3 

solver choice for velocity components jacobiuvw 
solver choice for pressure petsc 

parallel strategy for pressure corrector: matrixlevel 

initialize u field 1. 
initialize k field .1 
initialize e field .1 

Figure 9. Sample nphase.dat file 
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The second set of files are the grid files, written in COBALT unstructured format, 
cobalt.inp and cobalt.be. The cobalt. inp file completely defines the input grid in faced 
based fonnat. The file is written in COBALT format, as defined in Figure 10. The user 
need not concern himself with the content of this file although if a formatted cobalt.inp file 
is written by the grid generator, it will be human readable. 


ndim nzones nbc 

nvert nface ncell maxppf maxfpc 

xi yi zi 

x 2 y 2 z 2 


invert Ynyert invert 

nvfi (f_vi(ivf),ivf=l,nvfi) f_el 

f e2i 

nvf 2 (f_v 2 (ivf),ivf=l,nvf 2 ) f_el 2 

f_e2 2 

nvfnface (f_V n face(ivf),ivf= 1 ,nvf n fa C e) 

f_6 1 nface f_^2 n f ace 

Figure 10. cobalt.inp grid file format 


In this file the parameters are defined as follows: 

• ndim = # of dimensions always 3 for NPHASE-PSU, even if a 2D case is 

set up (see “Running Two-Dimensional Problems,’ below) 

• nzones = 1 (always for NPHASE-PSU) 

• nbc = total number of boundary conditions that are defined in cobalt.inp and 
cobalt.be 

• nvert = number of vertices in model 

• nface = number of vertices in model 

• ncell = number of vertices in model 

• maxppf = maximum number of vertices per face in any one face in domain 

• maxfpc = maximum number of faces per element in any one element in 
domain 

• Nvert yivert divert = vertex coordinates for vertex ivert 

• nvfiface = number of vertices on face iface 

• f_Viface = list of vertex numbers for face iface, listed in order around the 
periphery of the face such that the right-hand-rule applied to this ordering 
defines a direction from bounding element ) f_el to f_e2 
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• f_e 1 iface = element number on “low” side of face iface 

• f_e2if ace = element number on “high” side of face iface. If f_e2if ace < 0 then 
iface is a boundary face and f_e2if ace specifies the (negative of the) boundary 
identifier as defined in cobalt.be. 

The cobalt.be file defines boundary condition names. This ASCII file is written in 
standard COBALT boundary condition file format, an example of which is shown in Figure 
1 1 . The user need not concern himself with the content of this file unless it is not written by 
the grid generator (e.g., HARPOON), in which case the user must build this file by hand to 
confonn to the boundary condition numbering in the grid generator and thereby the 
cobalt. inp file. 

################################################## 

Boundary Condition Specification File for: 

Gridgen grid exported : Thu May 18 10:51:30 2006 

################################################## 

11 

WallOO 

Gridgen be region: 1 1 

Methods: User Created BC 

User data supplied here - see COBALT doc! 

################################################## 

9 

InflowOO 

Gridgen be region: 9 

Methods: User Created BC 

User data supplied here - see COBALT doc! 

################################################## 

10 

PressureOO 

Gridgen be region: 1 0 

Methods: User Created BC 

User data supplied here - see COBALT doc! 

################################################# 

Figure 11. cobalt.bc file format 


In this file three boundary conditions have been defined, Wall OO, Inflow OO and 
Pressure OO. This file tells the pre-processor, fump, that Wall OO faces in the cobalt. inp file 
will have identifier f_e2if ace = -11. Also, for lnflow_00, f_e2if ace = -9, and for Pressure_00, 
f_e2 iface -10. 

The third set of files is solution restart files. NPHASE-PSU always generates a 
restart file for each processor after execution completion (and optionally at intermediate 
iterations/timesteps as defined in “Control Commands” section below). These output files 
confonn to the naming convention nphase_restart_outxxx where xxx is the 3-digit processor 
identifier (range from 000 to 999). Solution restarts are invoked by: 1) copying each of 
these output restart files to a corresponding input restart file 

(i.e., mv nphase restart outOOO nphase restart inOOO) and then 2) activating the keyword 
“initialize run with restart file” in nphase.dat (i.e., uncomment by removing leading “#”) 
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Boundary Condition Specification 

NPHASE-PSU supports no-slip wall, porous wall, pressure, inflow, symmetry, 
far field and other specialized boundary conditions. The user needs to define boundary 
patch names that conform to what NPHASE-PSU is expecting. Specifically, NPHASE- 
PSU requires boundary patch names that conform to the following syntax: 
Boundarytypeboundarynumber, where Boundary type is either of these character strings: 
Wall, Porwall, Inflow, Symmetry, Farfield, Pressure. Boundarynumber is a 2-digit integer 
starting at 00. So for example there may be two inflows in a model with different attributes, 
these would be named by the user Inflow OO and InflowO 1 . These strings must appear in 
cobalt.be. The grid generators GRIDGEN, Pointwise and ICEM automatically propagate 
these names into cobalt.be upon grid output, provided the user names them in the grid 
generator. If HARPOON is used the user must define these in cobalt.be. 

Each boundary type has its own attributes, which are defined in nphase.dat as 
specified in the “Control Commands” section below. 

fump 

The pre-processor to NPHASE-PSU, fump, reads the cobalt. inp file (formatted or 
unformatted) and the cobalt.be files as input and performs two tasks: 

1) Executes domain decomposition by invoking METIS (2006) . fump does this by 
first extracting graph information for METIS. Specifically, each element in the 
grid is designated as a “vertex” in the graph and each element with which it 
shares a face represents an “edge”. The kmetis module is used to partition the 
graph into approximately equal sizes. 

2) Builds the internal pointer connectivity (e.g., facc^ clement, edge^ vertex, 
edge^ face), boundary pointer connectivity and interprocessor communication 
information needed by NPHASE-PSU 

fump is run interactively (or in a script) by simply typing the absolute path of the 
executable, fump, to a UNIX shell prompt, fump has four small user inputs: 1) # of 
processors to use for domain decomposition, 2) whether METIS is to be used or the user 
will be specifying a decomposition graph, 3) any scale factor the user may wish to apply to 
the grid, 4) whether or not a global vertex connectivity file is generated (specialized option 
- most applications set = 0). 

The output of fump is a series of files, unphase. gridxxx, where where xxx is the 3- 
digit processor identifier (range from 000 to 999). Each of these files is read by the 
corresponding processor from the executing front end of NPHASE-PSU. 

fump is written in ANSI-C and C++ and compiles (at least) under the gnu C 
compilers (gee, g++). The executable is generated by invoking make from the FUMP 
directory delivered with the software. The METIS libraries that fump requires are delivered 
with the software. 

Code Execution 

If a single processor job is required, one simply need type the name of the absolute 
path to the executable, nphase, into a UNIX shell prompt in the working directory. 
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NPHASE-PSU is instrumented with mpi for inter-processor communication. MPICH2 is 
used, so mpirun is executed behind the scenes by the invoked mpiexec wrapper and shown 
below. On most production UNIX cluster systems interactive invocation of mpirun and 
mpiexec are not allowed, rather, the user must build a submit script, and submit the job to a 
queueing system such as PBS. 

A typical run script for execution of NPHASE-PSU on the gale cluster at Penn State 
ARL is included in the tutorial section below. The input files (nphase.dat, unphase. gridxxx 
and nphase restart inxxx) must be located in the working directory from where the 
executable is invoked. Output files (discussed in next section) are written to this directory 
as well. 

Postprocessing 

Figure 12 illustrates the back end for NPHASE-PSU. 


nphaserestartoutxxx 


nphaserestartoutOOO 


301 


302 



ENSIGHT files (one 
set per processor) 


emerge or 
emergetrans 


merged ENSIGHT files 


Figure 12. Sketch of back end for NPHASE-PSU 
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As with the front end, the backend now also supports arbitrary polyhedra. 
Specifically, ENSIGHT GOLD fonnat is used to write the output files. ENSIGHT Versions 
8.0 and later will read and display these files. 


Output Files 


There are five classes of output files to NPHASE-PSU. The first is the standard 
ASCII printed output file, nphase.out, an example of which is included in Figure 13. This 
contains an echo of the input, residual history, and, if enabled in backend, printed field data 
for single processor jobs, (these print nodal commands are commented out in backend.c 
since it is unlikely that a user would ever wish to obtain a printed output for an unstructured 
domain). 
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NN NN NN 
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HH HH 
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AAAAAAAAA 
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EEEEEEEEE 

EEEEEEEEE 

EE 

EE 

EEEE 

EE 

EE 

EEEEEEEEE 

EEEEEEEEE 


NPHASE : A COMPUTER PROGRAM FOR THE PREDICTION OF MULTIFIELD 
FLOWS WITH MASS, MOMENTUM AND ENERGY TRANSFER 


Developed 
*** entered 
*** exiting 
*** entered 
*** exiting 
*** entered 


by: Rob Kunz 

pre_author: reading nphase.dat input *** 
pre_author: finished pass 1 on input file *** 
pre_author2: reading nphase.dat input *** 
pre_author2: finished pass 1 on input file *** 
author: reading nphase.dat input file *** 


iterations to perform 5 

number of fields 1 

time accurate simulation 

temporal discretization momentum 1 

physical timestep in seconds . 1 

number of physical timesteps 4 

transient file write frequency 1 

produce ensight output 

dont perform wall match logic 

inlet patch 1 0 

1.0 0 . 0 . 1 . 1 . 0.1 0.1 0 . 0 . 0 0 
pressure patch 1 0 
0 . 1 . 1 . 0.1 0.1 0 . 0 . 

turbulent flow high reynolds number k epsilon 

constant fluid molecular viscosity l.e-3 

constant fluid density 1 

function entry/exit echo off 

solver sweeps for u 3 

solver sweeps for v 3 

solver sweeps for w 3 

solver sweeps for p 6 

solver sweeps for k 3 

solver sweeps for e 3 

solver choice for velocity components jacobiuvw 
solver choice for pressure petsc 

parallel strategy for pressure corrector: matrixlevel 
initialize u field 1. 

initialize k field .1 

initialize e field . 1 


exiting author: finished 2nd pass on input file 


about to exit front end 


1 6 . 498e-02 0.000e+00 0.000e+00 1.515e+01 8.950e-03 5.938e-02 
1 1 . 168e-01 8 . 494e-02 9.137e-02 3.547e+00 2.788e-03 2.511e-02 
1 3 . 017e-02 2 . 215e-02 2.427e-02 1.142e+00 2.530e-03 1.011e-02 
1 1 . 281e-02 8 . 440e-03 9.386e-03 2.657e-01 1.725e-03 3.497e-03 
1 7 . 488e-03 3.947e-03 4.671e-03 1.383e-01 1.132e-03 1.649e-03 


-1 . 337e+01 
■3 . 696e+00 
■1 . 147e+00 
■4 . 204e-01 
■2 . 265e-01 


Figure 13. Sample nphase.out file 
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The second output file that is generated is standard output + standard error 
conventionally redirected to n.out, an example of which is included in Figure 13. This file 
contains a summary of the grid topology, front end progress, iteration history, user defined 
outputs, and summary flow rate and CPU performance information. 


12 processor run 
*** begin execution of nphase 


unphase. grid files generated by fump version 2.0 
nvert=409798 nnode=388114 nfedge=4332994 nface=1102458 


inlet . nbcfedge= 1024, inlet .nbcface= 256, 

partition . nbcfedge= 215690, partition . nbcface= 53192, 
pressure . nbcf edge= 2080, pressure .nbcface= 542, 
wall . nbcf edge= 116882, wall .nbcface= 32228, 

finished reading grid in read_grid_unstruct 


number of tets in model = 20732 
number of hexs in model = 320105 
number of prisms in model = 5442 
number of pyramids in model = 29983 
number of non-simplicial volume elements in model = 11852 


inlet . nbcf aceid= 1 

partition . nbcf aceid= 12 

pressure . nbcfaceid= 1 

wall . nbcf aceid= 5 


5. 0867325450934 e+ 04 


total mass of ifield = 0 in volume = 

solidflags: icountlO = 0 

Number of internal solid faces = 0 

about to exit front end 


time accurate simulation 

number of time steps to be run: 4 

number of pseudotimesteps per physical timestep: 

time step in seconds : 1 . 00000000e-01 


about to make subdirectory TRANSIN if it does not exist 
done with make subdirectory TRANSIN if it does not exist 


writing ensight output file(s) for current timestep 
finished writing ensight output file(s) for current timestep 


commencing inner iterations for t = 1 . 00000000e-01 

iphystimeglobal = 1 


iter fid ru rv rw 

1 1 6 . 4 98e-02 0.000e+00 0.000e+00 1. 

2 1 1 . 168e-01 8 . 4 94e-02 9.137e-02 3. 

3 1 3 . 017e-02 2.215e-02 2.427e-02 1. 

4 1 1 . 281e-02 8.440e-03 9.386e-03 2. 


515e+01 8 . 950e-03 5 
547e+00 2 . 788e-03 2 


142e+00 2. 
657e-01 1. 


530e-03 1 
725e-03 3 


re gmass 

. 938e-02 -1.337e+01 
. 511e-02 -3 . 696e+00 
. 011e-02 -1 . 147e+00 
. 4 97e-03 -4.204e-01 


ws/ cif 
. 659e-06 
. 435e-06 
. 455e-06 
. 453e-06 


1 7 . 488e-03 3.947e-03 4.671e-03 1.383e-01 1.132e-03 1.649e-03 -2.265e-01 3.455e-06 


writing ensight output file(s) for current timestep 
finished writing ensight output file(s) for current timestep 
total mass of ifield = 0 in volume = 5 . 0867325450934e+04 


commencing inner iterations for t = 
iphystimeglobal = 


2 . 00000000e-01 
2 




Figure 14. Sample n.out file (part 1) 
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commencing inner iterations for t = 4 . 00000000e-01 

iphystimeglobal = 4 


iter 

fid 

ru 

rv 

rw 

rp 

rk 

re 

gmass 

ws/ cif 

16 

1 

6 . 074e-03 

3 . 018e-03 

3 . 006e-03 

2 . 330e-01 

3 . 567e-03 

5 . 877e-03 

-2 . 64 6e-01 

3 . 485e-06 

17 

1 

1 . 721e-03 

6 . 843e-04 

6 . 456e-04 

1 . 321e-01 

1 . 260e-03 

3 . 092e-03 

-1 . 695e-01 

3 . 470e-06 

18 

1 

1 . 002e-03 

4 . 132e-04 

3 . 987e-04 

5 . 940e-02 

6 . 373e-04 

2 . 435e-03 

-6 . 422e-02 

3 . 438e-06 

19 

1 

5 . 080e-04 

2 . 381e-04 

2 . 078e-04 

2 . 989e-02 

3 . 783e-04 

1 . 976e-03 

-4 . 274e-02 

3 . 4 69e-06 

20 

1 

3 . 505e-04 

1. 604e-04 

1. 606e-04 

1 . 385e-02 

2 . 450e-04 

1 . 577e-03 

-1 . 894e-02 

3 . 429e-06 

total mass 

of ifield = 0 in 

volume = 

5.0867325450934e+04 







writing output restart file(s) at iter = 20 
finished writing output restart file(s) at iter = 20 


writing ensight output file(s) for current times tep 
finished writing ensight output file(s) for current timestep 

total wall secs elapsed = 3 . 00361910e+01 

total cpu secs used by all processors = 3 . 29250000e+02 

total elements on all processors = 388114 

cpu secs/element/iteration = 4 . 24166611e-05 

wall secs/element/iteration = 3 . 86950625e-06 

average wall secs/element for each iteration = 3 . 47251214e-06 

ratio of partition faces to elements = 1 . 37052521e-01 

about to enter back end 

writing ensight output file(s) 
finished writing ensight output file(s) 

Number of Inlet Boundaries = 1 

Area and Average Pressure for Inlet Boundary ID 0 = 1024 m**2, 0 Pa 
Mass Flow for Field 0 Through Inlet Boundary ID 0 = -1024 Kg/s 
Mass Flow for Field 0 Through All Inlet Boundaries = -1024 Kg/s 

Number of Pressure Boundaries = 1 

Area and Average Pressure for Pressure Boundary ID 0 = 1024 m**2, 0 Pa 
Mass Flow for Field 0 Through Pressure Boundary ID 0 = 1023.98 Kg/s 
Mass Flow for Field 0 Through All Pressure Boundaries = 1023.98 Kg/s 

Done with boundary mass flow print 

ifield, (massin-massout) /massin : 

*** ending execution of nphase ** 


0 1.8497181118102e-05 


Figure 15. Sample n.out file (continued) 


The third output file that is generated is resid.print which is text file containing only 
the residual history for convenient plotting. An example is included in Figure 15. Residuals 
in the code are not “true” residuals, rather they are the RMS of Acj), the change in value of 
the variable solved for at a given iteration. (Note that this is the value solved for by 
NPHASE-PSU, see equation (76) in the theory manual.). The columns in resid.print (and 
the residual prints in nphase. out and n.out for that matter) are column 1 : iteration number; 
column 2: field number (there will be a residual for each field [say liquid and gas]; columns 
3-6: RMS(Au), RMS(Aw), RMS(Aw), RMS(Ap); For multi-field simulations column 7 
contains the volume fraction residual for each field, RMS(Aa); For diabatic simulations 
column 8 contains the enthalpy residual, RMS(Ah); For two-equation turbulence 
simulations columns 9 and 10 contains the k and s residuals, RMS(Ak), RMS(As) (or q 
and co, k and co, for other turbulence models); column 1 1 contains the wall-clock seconds 
per (cell iteration* field) for that iteration. 


NASA/CR— 2014-2 1665 1 


41 

83 


1 

1 

1.996e-03 

2 

1 

1.968e-03 

3 

1 

1.416e-03 

4 

1 

1.095e-03 

5 

1 

8.819e-04 

G 

1 

7.288e-04 

7 

1 

6.596e-04 

8 

1 

G.206e-04 

9 

1 

5*890e-04 

10 

1 

5.605e-04 

11 

1 

5.347e-04 

12 

1 

4.976e-04 

13 

1 

4.812e-04 

14 

1 

4,669e-04 

15 

1 

4,548e-04 

16 

1 

4.404e-04 

17 

1 

4.198e-04 

1R 

1 

3 qqip-ri4 


0. 000e+00 0.000e+00 
1.736e-04 1.733e-04 
1.886e-04 1.883e-04 

1. G14e-04 1.610e-04 
1.336e-04 1.332e-04 
l,15Se-04 1.153e-04 
1.026e-04 1.023e-04 
9.439e-05 9.413e-05 
8.817e-05 8.796e-05 
8.313e-05 8.297e-05 
7.913e-05 7.901e-05 
7.523e-05 7.513e-05 
7.881e-05 7.868e-05 
7.254e-05 7.225e-05 
G,911e-05 G.875e-05 
6.653e-05 G.614e-05 
6.422e-05 S.381e-05 
k iq^o-m r iqip-rfi 


1.138e-02 0.000e+00 
9,678e-03 0.000e+00 
2.501e-02 0.000e+00 
2.79Ge-02 0.000e+00 
5.263e-02 0.000e+00 
l.G91e-02 0.000e+00 
1.255e-02 0.000e+00 
3.820e-03 0.000e+00 
1.047e-03 0.000e+00 
2.504e-03 0.000e+00 
1.099e-02 0.000e+00 
9.648e-02 0.000e+00 
G.25Ge-03 0.000e+00 
3.667e-03 0.000e+00 
1.714e-03 0.000e+00 
8.748e-04 0.000e+00 
3.219e-03 0.000e+00 
i laqo-n? n nniwnn 


0,000e+00 6.205e-03 
0.000e+00 5.294e-03 
0.000e+00 3.857e-03 
0.000e+00 3.151e-03 
0.000e+00 2.715e-03 
0.000e+00 2.403e-03 
0.000e+00 2.164e-03 
0.000e+00 1.977e-03 
0.000e+00 1.826e-03 
0.000e+00 1.703e-03 
0.000e+00 1.600e-03 
0.000e+00 1.512e-03 
0.000e+00 1.436e-03 
0.000e+00 1.370e-03 
0,000e+00 1.311e-03 
0.000e+00 1.259e-03 
0.000e+00 1.213e-03 
n nnn<.+nn i i7np-r>3 


2.659e-01 -7.525e-04 
5.72Ge-02 -1.426e-03 
4.036e-02 -1.723e-03 
3.965e-02 -1.837e-03 
3.441e-02 -1.635e-03 
2.881e-02 -1.195e-03 
2.407e-02 -7.996e-04 
2.025e-02 -2.927e-04 
1.719e-02 1.957e-04 
1.473e-02 6.738e-04 
1.275e-02 1.0Gle-03 
1.113e-02 7.315e-04 
9.803e-03 4.853e-04 
8.714e-03 2.734e-04 
7.812e-03 9.532e-05 
7.0G0e-03 -6.927e-05 
G.426e-03 -2.094e-04 
q qqq<.-riq -? Rq7o-n4 


Figure 16. Sample resid.print file 


The fourth set of output files are solution restart files. NPHASE-PSU always 
generates a restart file for each processor after execution completion (and optionally at 
intermediate iterations/timesteps as defined in “Control Commands” section below). These 
output files confonn to the naming convention nphase restart outxxx where xxx is the 3- 
digit processor identifier (range from 000 to 999). Solution restarts are invoked by 1) 
copying each of these output restart files to a corresponding input restart file (i.e., mv 
nphase restart outOOO nphase restart inOOO) and then 2) activating the keyword “initialize 
run with restart file” in nphase.dat. 

The fifth set of files are ENSIGHT output files. For steady state cases, each 
processor generates at least a geometry/grid file: engold.geo.xxx, pressure and velocity 
scalar files: engold.Esca.pOO.xxx, engold.Esca.un 77 .xxx, engold.Esca.vnn.xxx, and 
engold.Esca.wnn.xxx, where xxx is the processor number (0-999) and nn are the fields 
present (00-99). (Since NPFIASE-PSU employs a single pressure fonnulation only 
engold.Esca.pOO.xxx is generated). In addition to velocity and pressure files, other 
ENSIGHT files are generated depending on problem type and user specification (keyword 
- see Control Commands section below). Specifically, volume fraction files are generated 
for multi-field simulations: engold.Esca.ann.xxx; enthalpy and temperature files are 
generated for diabatic simulations: engold.Esca.hnn.xxx, engold.Esca.tnn.xxx ; turbulence 
quantity files are generated for turbulence simulations as appropriate for the model used: 
engold.Esca.knn.xxx (turbulent kinetic energy), engold.Esca.enn.xxx(turbulent dissipation 
rate, s), engold.Esca.mnn.xxx (eddy viscosity), engold.Esca.uunn.xxx, 
engold.Esca.vvnn.xxx, engold.Esca.wwnn.xxx, engold.Esca.uvnn.xxx, 
engold.Esca.vwnn.xxx, engold.Esca.wunn.xxx, (Reynolds stresses) for full Reynolds Stress 
models, engold.Esca.fnn.xxx (v2f damping function) and engold.Esca.vvnn.xxx (wall 
normal Reynolds stress) for the v2f model, and engold.Esca.nnn.xxx (Spalart-Allmaras 
kinematic-eddy-viscosity-like variable) for the Spalart-Allmaras model. The user can select 
two popular flow visualization parameters for ensight file generation: 
engold.Esca.helnn.xxx (relative helicity density) and engold.Esca.iswnn.xxx (intrinsic 
swirl). For compressible flow simulations, Mach number files are produced: 
engold.Esca.machnn.xxx. For simulations where interfacial area density transport is 
invoked, A’” (see eq. (7)) files are produced: engold.Esca.ainn.xxx. 
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emerge and emergetrans 


Once NPHASE-PSU has completed there exist many ENSIGHT files in the 
working directory from which the case is executed. One runs emerge to merge the 
ENSIGHT files generated on the different processors into single ENSIGHT files for each 
variable. One runs emergetrans to merge the ENSIGHT files generated on the different 
processors into single ENSIGHT files for each variable at each timestep that the user has 
selected for outputting these files, emerge and emergetrans are delivered in software 
directory EMERGE. They are written in FORTRAN and some C (I/O functions) and are 
compiled using makefiles (makeem, makeemt) that expect pgf90 and gcc, for which these 
codes are guaranteed to compile and run. Other compilers are untested. 

To run emerge, simply types “emerge” to the UNIX shell. The user will be 
prompted for the number of processors and the number of fields. Once emerge has 
completed, the following files are resident in the user’s working directory: engold.geo, 
engold.Esca.pOO, engold.Esca.unn, engold.Esca.v««, and engold.Esca.w«« (where nn are 
the fields present (00-99)). Other case specific merged files are also created (e.g., 
engold.Esca.aw7, engold.Esca.h 77 n, eng 0 ld.Esca.kw 7 , etc...). Finally emerge generates 
velocity vector files: engold.Evec.uvwnn. All files generated appear in the file engold.case, 
an example of which appears in Figure 17. Upon execution of emerge the user needs to 
transfer the case file, and all of the files it points to, to the file system where he will run 
ENSIGHT. To run ENSIGHT, the case file is read in and each of the other files are 
automatically brought in. 

# BOF: engold.case 

FORMAT 

type: ensight gold 

GEOMETRY 

model: engold.geo 

VARIABLE 


scalar 

per 

element : 

Pressure 


engold. pOO .Esca 

scalar 

per 

element : 

X-velocity 

O 

O 

engold. uOO .Esca 

scalar 

per 

element : 

Y-velocity 

O 

O 

engold. vOO .Esca 

scalar 

per 

element : 

Z-velocity 

O 

O 

engold . wOO . Esca 

scalar 

per 

element : 

TKE 

00 

engold. kOO .Esca 

scalar 

per 

element : 

TDS 

'00 

engold. eOO .Esca 

scalar 

per 

element : 

Eddy-viscosity 

O 

O 

engold. mOO .Esca 

vector 

per 

element : 

Velocity- vector 

_00 

engold . uvwOO . Evec 

# EOF: 

engold. case 





Figure 17. Example engold.case file 


To run emergetrans, simply type “emergetrans” to the UNIX shell. The user will be 
prompted for the number of processors, the number of fields, whether the grid is stationary 
or moving, the first integer timestep counter, the last integer timestep counter and the 
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interger timestep increment, emergetrans reads all of the transient per-processor files from 
the subdirectory, TRANSIN, and writes all output files to a subdirectory, TRANSOUT. 
Once emergetrans has completed, the following files are resident in TRANSOUT: 
engold.geo (or engold.geo.////// (if grid is moving), engold. Es ca.p 0 O.tttttt, 
engold.Esca.u nn.tttttt, engold.Esca ,\nn. tttttt, and engold. Esca. WW7. ////// (where nn are the 
fields present (00-99), and tttttt is the integer timestep). Other case specific merged files are 
also created (e.g., engold.Esca. ann. tttttt, engold.Esca. linn.//////, engold.Esca.knn.//////, 
etc...). Finally emergetrans generates velocity vector files: engold.Evec.uvw nn.tttttt. All 
files generated appear in the file TRANSOUT/engold transient.case, an example of which 
appears in Figure 18. Upon execution of emerge the user needs to transfer the case file and 
all of the files it points to the file system where he will run ENSIGHT. To run ENSIGHT, 
the case file is read in and each of the other files are automatically brought in. Upon 
execution of emergetrans the user needs to transfer the case file, and all of the files it points 
to, to the file system where he will run ENSIGHT. To run ENSIGHT, the case file is read 
in and each of the other files are automatically brought in. 


# BOF: 

engold transient . case 


FORMAT 




type: 

ensight 



GEOMETRY 



model : 

engold.geo 



VARIABLE 



scalar 

per element: 

Pressure 

engold. p00 .Esca. ****** 

scalar 

per element: 

X-velocity 00 

engold. uOO .Esca.****** 

scalar 

per element: 

Y-velocity 00 

engold.vOO. Esca. ****** 

scalar 

per element: 

Z-velocity 00 

engold.wOO. Esca. ****** 

vector 

per element: 

Velocity 00 

0 gold.uvwOOO vec . ****** 

scalar 

per element: 

TKE 00 

engoldOkOO .Esca.****** 

scalar 

per element: 

TDS 00 

engoldOeOO .Esca.****** 

scalar 

per element: 

Mu t_00 

engoldOmOO. Esca. ****** 

TIME 




time set: 1 



number 

of steps: 

5 


filename numbers: 




0 

1 2 

3 4 

time values: 



0 . 000E + 00 0 . 1 00E+01 0.200E + 01 

0 . 300E+01 0 . 400E+01 

# EOF: 

engold transient . case 




Figure 18. 

Example engoldtransient.case file 


Tutorials 

Four tutorials are provided as of Ver. 3.1, Rev. 1.12. The first is a general 
familiarization case, the second and third are specific to the DARPA FDR program which 
co-funded this documentation. The fourth is a highly compressible case. Each is delivered 
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Tutorials 

Four tutorials are provided as of Ver. 3.1, Rev. 1.12. The first is a general 
familiarization case, the second and third are specific to the DARPA FDR program which 
co-funded this documentation. The fourth is a highly compressible case. Each is delivered 
with the software and they unpack to ./NPHASE-PSU/TUTORIAL_l, ./NPHASE- 
PSU/TUTORIAL HIPLATE, ,/NPHASE-PSU/TUTORIAL_5415, ,/NPHASE- 
PSU/TUTORIAL_MACH_BUMP. Specifically, the necessary input files are provided for 
each tutorial: nphase.dat, cobalt. inp, cobalt.be, and a few other utility files. The user can 
follow the tutorials below successively applying the NPHASE-PSU pre-processing, 
execution and postprocessing steps, starting with these files. 

Tutorial Case 1: Turbulent, unsteady, arbitrary polyhedra two-body model 

This purely notional case is for flow around a deformed sphere in the vicinity of a 
solid wall in a channel. The grid was generated by HARPOON and is hex dominant, but 
contains “hanging nodes” which gives rise to elements that NPHASE-PSU interprets as 
arbitrary n-sided-polyhera (i.e., not tetrahedral, hexehedra, prisms or pyramids). The run is 
turbulent and unsteady. 

To start the tutorial the user needs to go to the TUTORIALl directory. There are 
six files there when the software is unpacked: nphase.dat, run.nphase, cobalt.inp, cobalt.be 
attributes. inp and mrf. The first step is to execute fump as illustrated in Figure 18. Here we 
have chosen 8 processors and no scaling of the geometry. This step generates the 8 files 
unphase. gridOOO, unphase.gridOOl, ...., unphase. grid007. Figure 19 shows the file 
nphase.dat. This is an unsteady run, with four timesteps. 5 inner iterations per timestep and 
an ENSIGHT output after every timestep. 

fump 

*** begin execution of fump *** 
enter number of processors for domain decomposition: 

12 

model will be decomposed into 12 partitions 
enter scale factor: 

1.0 

grid will be scaled by 1 
enter : 

0) metis decomposition : 

1) decomposition specified in graph.dat: 

0 

decomposition will be performed using metis 
in problemsize 
cobalt.inp is ascii 

*file_type 0 
ibc = 0, is a Wall OO boundary 
ibc = 1 , is a Inlet OO boundary 
ibc = 2, is a Pressure OO boundary 
ndim,nzones,nbc 3 13 

*nvert,*nface_cobalt,*ncell,maxppf,maxfpc = 40670 277884 125058 4 6 

*nface = 267148, pressure. nbcface = 1420, pressure. nbcfedge = 4260 
inlet.nbcface = 400, inlet.nbcfedge = 1600 

farfield.nbcface = 0, farfield.nbcfedge = 0 

Figure 19. fump run stream for TUTORIAL l. Bold-black=user supplied, blue=code generated. 
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#case title: 

#2 -CELL HANGING NODES 

iterations to perform 5 

number of fields 1 

time accurate simulation 
temporal discretization momentum 1 
physical timestep in seconds . 1 
number of physical timesteps 4 
transient file write frequency 1 

#initialize run with restart file 

produce ensight output 

dont perform wall match logic 

inlet patch 1 0 

1.0 0 . 0 . 1 . 1 . 0.1 0.1 0 . 0.0 0 

pressure patch 1 0 

0 . 1 . 1 . 0.1 0.1 0 . 0 . 

turbulent flow high reynolds number k epsilon 

constant fluid molecular viscosity l.e-3 
constant fluid density 1 

function entry/exit echo off 

solver sweeps for u 3 
solver sweeps for v 3 
solver sweeps for w 3 
solver sweeps for p 6 
solver sweeps for k 3 
solver sweeps for e 3 

solver choice for velocity components jacobiuvw 

solver choice for pressure petsc 

parallel strategy for pressure corrector: matrixlevel 

initialize u field 1 . 
initialize k field . 1 
initialize e field . 1 

Figure 20. nphase.dat file for TUTORIAL_l. 


The file run.nphase is the job submit script that needs to be edited by the user to 
confonn to the system he is running on. The delivered version appears in Figure 20, and is 
that used to execute on Penn State ARL clusters UWE, STU and SADIE. The user must 
edit the script to point to mpiexec and the executable (highlighted in blue) on their system, 
as well as the nodes=x:ppn=x line to conform to the number of nodes and processors for 
the current run. Standard PBS attributes are included that presumably are machine 
independent. 

The job is submitted to PBS using the qsub command: qsub run.nphase (unless 
there are other submit scripts available on the user’s system). This 388,1 14 element 12 
processor job runs in about 90 seconds on a modern LINUX cluster. The code generates 
many files, most of them ENSIGFIT files, as can be verified by doing an “Is”. 
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#!/bin/bash 

#PBS -1 nodes=l :ppn=12 

#PBS -1 walltime=l :00:00 
#PBS -j oe 
#PBS -q batch 

cd $PBS O WORKDIR 

export PATH=/home/local/mpich2-l ,0.5/intel/bin:$PATH 

# Calculate the number of processors 
NPROCS='wc -1 < $PBS NODEFILE' 
echo SNPROCS 

# Create nodefile (with duplicates removed) for use by mpdboot 
NODEFILE_UNIQ=/tmp/'basename $ {PBS NODEFILE}' .uniq 
NODEFILE _MPDS=/tmp/'basename $ {PBSNODEFILE} ' .rmpd 
cat SPBSNODEFILE j uniq > SNODEFILE UNIQ 
N_MPDS='wc -1 < SNODEFILE UNIQ ' 

tail -"$(($N_MPDS- 1 ))" SNODEFILEUNIQ > SNODEFILE MPDS 

echo "all nodes:" 
cat SNODEFILE UNIQ 
echo "remote nodes:" 
cat SNODEFILE MPDS 

# Boot the MPI2 engine. 

/home/local/mpich2-1.0.5/intel/bin/mpdboot — totalnum=$N_MPDS -verbose — rsh=ssh — file=${NODEFILE_MPDS} 
#/home/local/mpich2- 1 ,0.5/intel/bin/mpdtrace 

# Run your executable 


The user can view nphase.out, n.out (which is the standard output and standard 
error redirect of the NPHASE-PSU run,) and resid.print. The next postprocessing step is to 
execute emerge which produces processor merged ENSIGHT output files of the final 
solution. The emerge run stream for this case is included in Figure 21. Once emerge is 
completed the user can migrate the merged data files and case files to the location he plans 
on postprocessing using ENSIGHT. These files are: engold.case, engold.geo, 
engold.wOO.Esca, engold.vOO.Esca, engold.uvwOO.Evec, engold.uOO.Esca, 
engold.pOO.Esca, engold.mOO.Esca, engold.kOO.Esca, engold.eOO.Esca. The user can also 
execute emergetrans for this case which produces processor merged ENSIGHT output files 
at every timestep. The emergetrans run stream for this case is included in Figure 22. Once 
emergetrans is completed the user can migrate the merged data files and case files to the 
location he plans on postprocessing using ENSIGHT. These files are the complete contents 
of the subdirectory TRANSOUT. 



/home/local/mpich2-1.0.5/intel/bin/mpiexec -np SNPROCS ,/nphase > n.out 
/home/local/mpich2- 1 .0.5 int el/bin/mp dallexit 


Figure 21. run.nphase script for TUTORIAL l. 


NASA/CR— 2014-2 1665 1 


47 

89 


emerge 


egoldmerge assembly 

enter number of partitions/processors 

12 

enter number of fields 

1 

engold.geo.000 

1 49860 

1 numcellstets(iproc) = 2555 

1 numcellshexs(iproc) = 39823 

1 numcells_prisms(iproc) = 650 


no WU file 
reading F-V2F files 
no F_V2F file 
reading IB LANK files 
no IB LANK file 
reading Temperature files 
no Temperature file 
reading Aint files 
no Aint 

re-reading velocity files 
FORTRAN STOP 


Figure 22. emerge run stream for TUTORIAL!. Bold-black=user supplied, blue=code generated. 
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emergetrans 


egoldmerge assembly for transients 


enter number of partitions/processors 

12 

enter number of fields 

1 

enter: 

1) stationary mesh 

2) moving mesh 

1 

enter the integer number of the first timestep 

0 

enter the integer number of the last timestep 

4 

enter the integer timestep increment 

1 

enter: 

1) use nphase.out to decide which variables to process 

2) use efiles.dat to decide which variables to process 

1 

reading from nphase.out to determine which ensight 
transient files were written by run being postprocessed 
Ensight transient variable files to be merged 
Pressure 
X-velocity 


no VW file 
reading WU files 
no WU file 
reading F-V2F files 
no F file 

reading IB LANK files 
no I file 

reading Temperature files 
no Temperature file 
reading Aint files 
no Aint file 

re-reading velocity files 
FORTRAN STOP 


Figure 23. emergetrans run stream for TUTORIALl. Bold-black=user supplied, blue=code generated. 


The second part of this tutorial involves modifying nphase.dat to run a steady state 
problem, and then running this problem to convergence, including doing a mid-run restart. 
To do this copy nphase.dat to nphase.dat. sav, so it can be recovered later. Next change the 
number of iterations to 250 and comment out the unsteady run keywords as shown (in blue) 
in Figure 23. Submit and run the job to completion of 250 iterations. To execute a restart 
comment in the “initialize run with restart file” keyword (as shown in green in Figure 23), 
and copy all of the nphaserestartoutxxx files to nphase restart inxxx. A convenient way 
of doing this is to build a small script, mrf, as included in the TUTORIAL l directory and 
as appears in Figure 24. Resubmit, running the code out to 500 iterations, and then run 
emerge again. 
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Figure 26 and 27 show the convergence history (TECPLOT used with resid.print as 
input) and a representative ENSIGHT visualization of the output of this steady state run 


case title: 

2-CELL HANGING NODES 

iterations to perform 250 

number of fields 1 

#time accurate simulation 
#temporal discretization momentum 1 
#physical timestep in seconds . 1 
#number of physical timesteps 4 
transient file write frequency 1 

initialize run with restart file 

produce ensight output 


Figure 24. modified nphase.dat file for TUTORIALl. 


mv nphase restart outOOO nphase restart inOOO 
mv nphase restart outOOl nphase restart inOOl 
mv nphase_restart_out002 nphase_restart_in002 
mv nphase_restart_out003 nphase_restart_in003 
mv nphase_restart_out004 nphase_restart_in004 
mv nphase_restart_out005 nphase_restart_in005 
mv nphase_restart_out006 nphase_restart_in006 
mv nphase_restart_out007 nphase_restart_in007 

Figure 25. mrf script 



Figure 26. TECPLOT plot of residual history of steady TUTORIAL_l simulation. 
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!sBf 






s® s -«a 


Eddy_viscosityJX) 


4.748e-002 


3.56 e-002 


2.374e-002 

l.J87e-002 

O.OOOe+OOO 


X_velociiy_00 


1 .492e+0()0 
1.080e+000 
6.677e-00] 

2.556e-00J 


-i.565e-001 

Pressure 


5.303e-001 


2.132e-001 


-1.039e-001 

-4.210e-001 


-7.381e-001 


Figure 27. ENSIGHT visualization of steady TUTORIAL l simulation. 
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Tutorial Case 2: Multiphase HIPLATE Simulation 

This case demonstrates the application of NPHASE-PSU to the HIPLATE 
configuration. The HIPLATE experimental programs were carried out by Ceccio and his 
colleagues at the Uni vers ty of Michigan ( Sanders et ah, 2006) under the DARPA Friction 
Drag Reduction program. Comparison of NPHASE-PSU to and calibration of NPHASE- 
PSU against these data sets were elements of the program (Kunz, et ah, 2006 , Kunz et ah, 
2007) . HIPLATE is a very high Reynolds number plate configuration tested in the US 
Navy’s Large Cavitation Channel, using microbubble and polymer drag reduction schemes. 

Although NPHASE-PSU was validated against and/or applied to the entire 
spectrum of HIPLATE runs (five gas injection rates, three tunnel speeds, various injector 
port gas flow splits, with vs. without surfactant, two injector geometries, and two water 
tu nn el test programs [HIPLATE I, II]), we here demonstrate a coarse grid (fast running) 
application corresponding to the January 2005 program deliverable (validation of code 
against HIPLATE I data). Specifically, a 12 m/s, Q ga s=800 scfm case is set up and run. The 
physical model set used here is very close to that used for the January 2006 deliverable to 
DARPA which was a validation study comparing NPHASE-PSU simulations to all 
HIPLATE I data. 

To start the tutorial the user needs to go to the TUTORIAL HIPLATE directory. 
There are four files there when the software is unpacked: nphase.dat, run.nphase, 
cobalt. inp, and cobalt.be. The first step is to execute fump (as described in the 
Preprocessing section and Tutorial l section above), using (here) 1 domain and a scaling 
factor of 1.0. This step generates a single grid+topology file, unphase. gridOOO. Figure 28 
shows the file nphase.dat file used for this tutorial. The keywords “employ hiplate 
modeling” are included to instruct nphase to execute the function hiplate output.c each 
iteration which prints skin friction data (and, if commented in, other data) to standard 
output (redirected to n.out in run.nphase). (See the section entitled Building User 
Specific/Case Specific Postprocessing for details on how to modify or adapt this kind of 
output for different output or simulation case.) 


First the code has to be run in single phase mode to generate the comparison flat 
plate skin friction values at the 6 axial locations on the HIPLATE where shear stress 
measurements were made. To do this, the user needs to modify the green highlighted 
sections in Figure 28 as follows: 1) change number of iterations to 500, 2) change porous 
wall injection velocity from 9.0476 to 0.0000 and the porous wall permeability from 0.7 to 
0.0., 3) comment out the relaxation factor lines (allow for defaults). So the code is being 
run in full 2-phase mode but there is no gas present or injected to it returns a single phase 
convergence and results. 

This single phase job can be run interactively (since its only one processor) or 
submitted to PBS using the qsub command: qsub run.nphase (unless there are other submit 
scripts available on the user’s system). This single phase job processor job executes 500 
iterations in about 55 wall seconds on a modern LINUX cluster. Figure 29 shows a section 
of the n.out file for this single phase run that includes the computed Cf distribution at the 
six axial HIPLATE stations. To run the microbubble drag reduction case, the nphase.dat 
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iterations to perform 20000 

number of fields 2 

produce ensight output 

restart file write frequency 1000 
dont perform wall match logic 

employ hiplate modeling 

gravity vector 0. 9.81 0. 
employ modified pressure 

reference density for modified pressure 1 

overwrite inlet patch boundary conditions on restart 
overwrite porous wall patch boundary conditions on restart 
overwrite pressure patch boundary conditions on restart 
overwrite_gas_molecular_viscosities 200 1 . 5e+3 

inlet patch 2 0 

12.0 0. 0. 1000. 0.999999 1.33e-7 1.76e-6 0. 1. 

12.0 0 . 0 . 1 . 0.000001 0 . 0 . 0 . 1 . 


porous wall patch 2 0 


9.0476 

0. 

0. 1000. 

0.000001 1 . 33e-7 1.76e-6 0. 

0.7 

9.0476 

0. 

0. 1. 

0.999999 0. 0. 0. 

0.7 


pressure patch 2 0 

0.0 1000. 0.999999 1.33e-7 1.76e-6 0. 

0.0 1 . 0.000001 0 . 0 . 0 . 

constant fluid density 1000. 1. 

constant fluid molecular viscosity l.e-3 1.5e-5 
constant fluid surface tension .072 0. 

Turbulence model for each field 1 0 

function entry/exit echo off 
adiabatic flow 

solver sweeps for u 3 
solver sweeps for v 3 
solver sweeps for a 3 
solver sweeps for p 10 
solver sweeps for tke 3 
solver sweeps for tds 3 

solver choice for velocity components jacobiuvw 
solver choice for pressure petsc 

parallel strategy for pressure corrector: matrixlevel 
employ rhie chow for dispersion terms 


relaxation 

factor 

for 

u 

.05 

.05 

relaxation 

factor 

for 

V 

.05 

.05 

relaxation 

factor 

for 

a 

.05 

.05 


initialize 

initialize 

initialize 

initialize 

initialize 

initialize 


u field 12. 12. 
v field 0. 0. 
w field 0. 0. 
a field .999999 0.000001 
tke field 1.33e-7 0. 
tds field 1.76e-6 0. 


i f 
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interfacial area transport model for each field 0 1 
interfacial area coalescence model for each field 0 3 

initialize interfacial area using volume fraction and characteristic length 
interfacial area feeds back into bubble diameter frequency 100 
constant field characteristic diameter 99. .000400 
interfield drag models 1 
1 2 6 .8 

bubble cluster drag modification 
interfield nondrag models 5 
1 2 2 5 50. 

1 2 2 7 1. 2. -2.0 0.5 2. 

1 2 3 3 .1 
1 2 1 6 .25 
12531 . 


Figure 28. nphase.dat file for HIPLATE tutorial 


iter fid ru rv rw 

1 1 3 . 580e-03 4.631e-01 0.000e+00 

1 2 1 . 005e+00 1 . 672e-14 0.000e+00 

cfl, cf2, cf3, cf4, cf5, cf6 2 . 57 927225e-01 
clipcountO 64 

2 1 7 . 4 90e-02 9.948e-02 0.000e+00 

2 2 5 . 775e-01 9.039e-02 0.000e+00 

cfl, cf2, cf3, cf4, cf5, cf6 4 . 48393267e+00 


rp ra raint rh rk . . . 

5 . 767e+03 2.210e-08 0.000e+00 0.000e+00 1.104e-05 ... 
5 . 767e+03 2.139e-08 3.208e-04 0.000e+00 0.000e+00 ... 
2 . 57927225e-01 2 . 57927225e-01 2 . 57927225e-01 

1.157e+03 3 . 543e-08 0.000e+00 0.000e+00 1.324e-03 ... 
1.157e+03 3 . 467e-08 5.207e-04 0.000e+00 0.000e+00 ... 
5 . 93114814e+00 6 . 51944986e+00 6 . 64137946e+00 


4TO 1 8 . 279e-06 1.624e-07 0.000e+00 

499 2 5 . 797e-01 4.853e-01 0.000e+00 
cfl, cf2, cf3, cf4, cf5, cf6 1 . 5924 67 92e+02 

500 1 8 . 102e-06 1.540e-07 0.000e+00 

500 2 1 . 103e-01 1.703e-01 0.000e+00 

cfl,cf2,cf3,cf4,cf5,cf6 1 . 5924 67 94e+02 
final cfl,cf2,cf3,cf4,cf5,cf6 1 . 59246794e+02 
avgcf cfl,cf2,cf3,cf4,cf5,cf6 1 . 56126868e+02 


3 . 534e-02 6.881e- 
3 . 534e-02 4.581e- 
1 . 4 9083851e+02 
3 . 885e-02 3.922e- 


06 0 . 000e+00 0.000e+00 2.490e-06 
06 6 . 947e-02 0.000e+00 0.000e+00 
1 ,39654688e+02 1 . 36211604e+02 

06 0 . 000e+00 0.000e+00 2.440e-06 


3 . 885e-02 3.541e-06 5.406e-02 0.000e+00 0.000e+00 
1 . 49083857e+02 1 . 39654745e+02 1 . 36211706e+02 

1 . 49083857e+02 1 . 39654745e+02 

1 ,41946021e+02 1 . 39721464e+02 


stdv cfl, cf2, cf3, cf4, cf5, cf6 1 . 36051236e+01 1 . 79146591e+01 1 . 66518929e+01 


Figure 29. Snippets from standard output (n.out) for single phase TUTORIAL HIPLATE simulation, 
illustrating case specific skin friction coefficient output. 
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ru 


rv 


rw 


ra 


1 


iter fid ru rv rw rp ra raint rh rk 

1 1 2 . 007e-04 1.944e-02 O.OOOe+OO 7.098e+05 2.068e-03 0.000e+00 O.OOOe+OO 1.102e-05 

2 7 . 272e-02 1.021e-02 O.OOOe+OO 7.098e+05 2.067e-03 9.201e+03 O.OOOe+OO O.OOOe+OO 


cfl,cf2,cf3,cf4,cf5,cf6 2 . 57927225e-01 2 . 57927225e-01 2 . 57927225e-01 2 . 57927225e-01 

clipcountO 85 

2 1 4 . 648e-01 2.057e-01 O.OOOe+OO 6.805e+05 1.970e-03 O.OOOe+OO 0, 

2 2 7 . 036e-02 3.113e-03 O.OOOe+OO 6.805e+05 1.969e-03 7.497e+03 0. 
cfl,cf2,cf3,cf4,cf5,cf6 4 . 53209563e+00 5 . 96370783e+00 6 . 54783080e+00 


000e+00 1 . 373e-03 
000e+00 O.OOOe+OO 
6 . 66829586e+00 


19999 1 1 . 4 63e-05 4.050e-06 

19999 2 1 . 963e-03 1.345e-04 

cfl,cf2,cf3,cf4,cf5,cf6 7.61 

20000 1 1 . 447e-05 4.117e-06 

20000 2 2 . 261e-03 1.369e-04 

cfl,cf2,cf3,cf4,cf5,cf6 7.61 
final cf 1 , cf2 , cf3 , cf 4 , cf5 , cf 6 
avgcf cf 1 , cf2 , cf3 , cf 4 , cf5 , cf 6 
stdv cf 1 , cf2 , cf3 , cf 4 , cf5 , cf 6 


O.OOOe+OO 8 . 241e-01 3.974e-04 0 
O.OOOe+OO 8 . 241e-01 3.906e-04 8 
6747 67e+01 9 . 95633653e+01 1.12 

O.OOOe+OO 8 . 499e-01 3.973e-04 0 
O.OOOe+OO 8 . 499e-01 3.905e-04 8 
674747e+01 9 . 95633580e+01 1.12 

7 . 61674747e+01 9 . 95633580e+01 

7 . 6167 6731e+01 9 . 95634673e+01 

1 . 64306879e-04 2 . 11569909e-04 


. 000e+00 O.OOOe+OO 1.243e-05 
. 455e+01 O.OOOe+OO O.OOOe+OO 
126232e+02 1 . 20645034e+02 

. 000e+00 O.OOOe+OO 1.226e-05 
. 445e+01 O.OOOe+OO O.OOOe+OO 
101175e+02 1 . 20644965e+02 

1 . 12101175e+02 
1 . 12181731e+02 
5 . 42004104e-02 


Figure 30. Snippet from standard output (n.out) for two-phase TUTORIAL HIPLATE simulation, 
illustrating case specific skin friction coefficient output. 
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Figure 31. Predicted skin friction coefficient vs. iteration at the six HIPLATE stations for the two- 

phase TUTORIAL_HIPLATE simulation. 
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Figure 32. Predicted skin friction coefficient vs. x at the six FilPLATE stations for the two-phase 

TUTORIAL HIPLATE simulation. 
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Figure 30 shows a section of the n.out file for the two-phase microbubble drag 
reduction run that includes the computed Cf distribution at the six axial HIPLATE stations. 

Figure 31 shows the computed 2-phase skin friction iteration history. This solution 
converged to the final Cf values in about 5000 iterations. Figure 32 shows the predicted 
drag reduction vs. x for this case. These results are nearly identical to the January 2006 
submission (grid is coarser here). Figure 33 shows contours of predicted gas volume 
fraction with the vertical axis scaled by a factor of 100 in order that the thin microbubble 
layer can be viewed. 


V olume_fraction_0 1 

l.OOOe+OOO 

7.500e-001 

5-000e-001 

2.500e-001 

2.603e-034 


Figure 33. Predicted contours of gas volume fraction for the two-phase TUTORIALHIPLATE 
simulation. Vertical axis scaled by a factor of 100. 
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Tutorial Case 3: Multiphase 5415 Simulation 

This case demonstrates the application of NPHASE-PSU to a Navy relevant 
configuration in 1 -phase and 2-phase modes. The geometry is the 5415 model, a hull form 
representative of the Arleigh-Burke class destroyer. This model has been extensively tested 
at the Naval Surface Warfare Center Carderock Division (http://www.dt.navy.mil/hyd/sur- 
shi-mod/index.html , for example.) The grid was generated usin g ICEM-CFD and is a 
standard tetra+prism mesh. Specifically, the mesh contains 1438852 elements of which 
463837 are prisms extruded in (nominally) five layers from the triangulated hull surface. 
This comparatively coarse mesh is suitable for wall function turbulence modelling. The 
hull surface incorporates a standard solid wall patch and a porous wall patch, aft of the bow 
dome, from which gas is injected. Figure 34 shows two views of the mesh used. 

To start the tutorial the user needs to go to the TUTORIAL 5415 directory. There 
are five files there when the software is unpacked: nphase.dat, run.nphase, cobalt. inp, 
cobalt.be and mrf. The first step is to execute fump (as described in the Preprocessing 
section and Tutorial _1 section above), using (here) 32 domains and a scaling factor of .00 1 
(to convert the grid from mm to m.) This step generates the 32 files unphase. gridOOO, 
unphase.gridOOl, ...., unphase.grid0031. Figure 35 shows the file nphase.dat. The 
keywords “employ model 5415 modeling” are included to instruct nphase to execute the 
function model5415_output.c each iteration which generates a wetted area and a net vehicle 
drag coefficient printout to standard output (redirected to n.out in run.nphase). (See the 
section entitled Building User Specific/Case Specific Postprocessing for details on how to 
modify or adapt this kind of output for different output or simulation case.) A freestream 
velocity of 2.2134 m/s is specified (model scale) and the k-s model is selected. 



Figure 34. Views of mesh employed for 5415 tutorial 


First the user should run the code to generate a single phase solution for comparison 
to multiphase runs to be made later. To do this, the user needs to modify the green 
highlighted sections in Figure 35 as follows: 1) change number of iterations to 2500, 2) 
change porous wall injection velocity from 0.3000 to 0.0000 and the porous wall 
permeability from 0.7 to 0.0. So the code is being run in full 2-phase mode but there is no 
gas present. The job is submitted to PBS using the qsub command: qsub run.nphase (unless 
there are other submit scripts available on the user’s system). This 32 processor job 
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executes 2500 iterations in about 20000 wall seconds on a modern LINUX cluster 
(banyan.dt.navy.mil). As before, the user can view nphase.out, n.out, and resid.print. Then 
emerge is run and the user can migrate the merged data files and case files to the location 
he plans on postprocessing using ENSIGHT. These files are: engold.case, engold.geo, 
engold.wOO.Esca, engold.vOO.Esca, engold.uvw00.Evec, engold.uOO.Esca, 
engold.p00.Esca, engold.mOO.Esca, engold.kOO.Esca, engold.eOO.Esca. 


#case title: 

#5415 

iterations to perform 10000 

employ model 5415 modeling 

number of fields 2 

#initialize run with restart file 

produce ensight output 

restart file write frequency 100 

dont perform wall match logic 
read wall proximity from file 


gravity vector 0. 0. -9.81 
employ modified pressure 

reference density for modified pressure 1000. 

#simple hydrostatic pressure treatment 
overwrite inlet patch boundary conditions on restart 
overwrite porous wall patch boundary conditions on restart 
overwrite pressure patch boundary conditions on restart 
#overwrite_gas_molecular_viscosities 200 1 . 5e+3 


inlet patch 2 0 

2.2134 0. 0. 1000. .99999999999 .002939 .02618 0. 0. 
2.2134 0. 0. 1. .00000000001 .0 .0 0. 0. 


pressure patch 2 0 

0.0 1000. .99999999999 .002939 .02618 0. 0. 

0.0 1 . .00000000001 .0 .0 0 . 0 . 


porous wall patch 2 0 


0.3000 

0. 0. 1000. 

.00000000001 

.002939 

.02618 0. 

0.7 

0.3000 

0. 0. 1. 

.99999999999 

.0 

o 

o 

0.7 


turbulence model for each field 1 0 


constant fluid molecular viscosity l.e-3 1.5e-5 
constant fluid density 1000. 1. 
constant fluid surface tension .072 0. 


spatial discretization momentum 2 


function entry/exit echo off 
adiabatic flow 

solversweepsforu 3 3 
solversweepsforv 3 3 
solversweepsforw 3 3 
solversweepsfora 3 3 
solversweepsfork 3 3 
solversweepsfore 3 3 
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solver choice for velocity components jacobiuvw 
solver choice for pressure petsc 

parallel strategy for pressure corrector: matrixlevel 
petscprintnorm 

employ rhie chow for dispersion terms 


#relaxation factor 
#relaxation factor 
#relaxation factor 
#relaxation factor 
#relaxation factor 
#relaxation factor 
# 

initialize u field 
initialize v field 
initialize w field 
initialize a field 
initialize k field 
initialize e field 


for u .2 .2 
for v .2 .2 
for w .2 .2 
for a .2 .2 
for k .2 .2 
for e .2 .2 

2.2134 2.2134 

0 . 0 . 

0 . 0 . 

.99999999999 .00000000001 
.002939 .0 
.02618 .0 


clips_on_velocity -5. 10. -5. 5. -5. 5. 

#interfacial area transport model for each field 0 1 
#interfacial area coalescence model for each field 0 3 

#initialize interfacial area using volume fraction and characteristic length 
#interfacial area feeds back into bubble diameter frequency 100 
constant field characteristic diameter 99. .000400 
interfield drag models 1 
1 2 6 0.8 

bubble cluster drag modification 
interfield nondrag models 5 
1 2 2 5 50. 

1 2 2 7 1. 2. -2.0 0.5 2. 

1 2 3 3 .1 
1 2 1 6 .25 
12531 . 


Figure 35. nphase.dat file for TUTORIAL_5415. 


Figure 35 shows the convergence history for this single phase case. The code 
converges well and then eventually flatlines due to weak unsteadiness in the solution. 

2500 1 4 . 207e-06 2.098e-06 8.752e-07 1.029e-03 2.231e-ll 0.000e+00 1.127e-06 2.070e-04 

2500 2 1 . 480e-01 2.008e-01 6.633e-02 1.029e-03 2.395e-ll 0.000e+00 0.000e+00 0.000e+00 

Awetted dot i: 2 . 2368721880023e-05 m A 2 

Viscous drag force: 3 . 3070651056143e+01 

Pressure drag force: 1 . 0161932291489e+01 

Total drag force: 4 . 3232583347 632e+01 

Figure 36 shows a section of the n.out file which illustrates the computed frontal 
wetted area and drag force computation. 

The second part of this tutorial involves running a notional microbubble drag 
reduction case, for comparison to the single phase drag just computed. To execute this part 
of the tutorial, the user should employ the nphase.dat file as it occurs in the software 
deliverable directory, shown in Figure 34. Numerous multifield extensions to a single 
phase run are apparent. Specifically: 

1) a wall proximity computation (for all cells not just wall adjacent cells) is 
automatically carried out to support some of the interfacial force models (wall-lift). This 
wall proximity computation can be quite time consuming the first run of this case (several 
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minutes). So the user should not be concerned. After a few minutes this process will be 
complete and the files wall_proximity.datax (xxx = processor number, 000-03 1 here) will 
appear in the working directory and the iterations will begin. For all subsequent runs using 
this grid the user can forego the time consuming wall proximity calculation by reading 
these files on input, by commenting in the keyword: “read wall proximity from file”. 

2) Gravity of course needs to be specified. Here a conventional modified pressure 
approach is used with a reference density set equal to the liquid density. This zeroes the 
magnitude of the buoyancy term in the liquid momentum equation which has been found to 
improve convergence. A simple hydrostatic head initialization and exit pressure boundary 
specification is invoked. 



Figure 36. TECPLOT plot of residual history of single phase TUTORIAL 5145 simulation. 


iter fid ru rv rw rp ra rh rk re 


2500 1 4 . 207e-06 2.098e-06 8.752e-07 1.029e-03 2.231e-ll 0.000e+00 1.127e-06 2.070e-04 ... 

2500 2 1 . 480e-01 2.008e-01 6.633e-02 1.029e-03 2.395e-ll 0.000e+00 0.000e+00 0.000e+00 ... 
Awetted dot i: 2 . 2368721880023e-05 m A 2 

Viscous drag force: 3 . 3070651056143e+01 

Pressure drag force: 1 . 0161932291489e+01 

Total drag force: 4 . 3232583347 632e+01 

Figure 37. Snippets from standard output (n.out) for single phase TUTORIAL5145 simulation, illustrating 

case specific drag output. 
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3) Inlet, pressure and porous boundary conditions have attributes specified for both 
fields. Note here that we specify liquid and gas volume fractions of lxlO' 6 within 0 and 1, 
again arising for robustness (matrix singularity) reasons. If a porous wall normal velocity 
of 0.38 m/s were set, this would yield a nondimensional gas injection rate, Q a /(Q a +Q w ) = 
0.44 which corresponds to the HIPLATE 12m/s, 800 scfm case analyzed in 
TUTORIALHIPLATE (based on a representative Q w =Ub(S-8 ), where 8 and 8* are 
computed based on standard flat plate momentum integral relations at the axial injector 
location). Accordingly a range of injection velocities ranging from 0.0 to 0.40 m/s were 
run, thereby spanning the Q a /(Q a +Q w ) range of all bubbly flow HIPLATE runs. The 
particular case detailed here is for Ui n j=0.3 m/s. 

4) Turbulence scalars are solved only for the continuous field here. 

5) Density and viscosity are set for both fields 

6) Solver sweeps are set for both fields for all scalars (except pressure). Relaxation 
factors are set for both fields for all scalars (except pressure). 

7) Initial values are set for both fields for all scalars. Initially the domain is filled 
with water (a =1, a~=0.) 

8) The remainder of the keywords specify various interfacial dynamics, mass 
transfer and interfacial area density transport models. 

The grid supplied in this tutorial is fairly coarse for a high Reynolds number 
boundary layer flow. Only five prism layers are employed and wall function resolution is 
specified. Accordingly the first order numerics associated with the volume fraction in this 
case gives rise to some numerical smearing of the gas layer. 

The model set chosen for this simulation includes all relevant interfacial dynamics 
(drag and non-drag) forces (see HIPLATE tutorial above and Kunz et ah, (2006)) . The 
interfacial force model set is identical to that presented for the HIPLATE tutorial. 

Typically, for microbubble drag reduction applications one runs the code until it 
“flatlines”. At that point one can look at the Cd history with iteration and its standard 
deviation. If the standard deviation is small compared to the magnitude of the predicted 
drag, then one accepts the mean Cd value as the prediction. 

In this Uj n j=0.3 simulation, the code was run for 10000 iterations. Figure 37 shows 
the first 2000 iterations of the convergence history for this two-phase case. The code 
flatlines due to gas velocity cutoffs that are specified (see Figure 34) which keep the gas 
velocities reasonable in regions where the gas fraction is vanishingly small (0(1 O' 6 )). 

Figure 38 shows the drag history for this case illustrating that for MBDR applications, 2500 
iterations is adequate for injection scheme assessment (and this is what was used for the 
other five 2-phase cases run as discussed below). Figure 39 shows a section of the n.out file 
which illustrates the computed frontal wetted area and drag force computation. Figure 40 
shows a view of the hull surface pressure and a gas volume fraction isosurface, a gas =0.5, 
illustrating the predicted topology of the bubble layer. 

The last set of results included in this tutorial parameterize the gas flow rate as 
might be done in a design assessment. Specifically, the code was run an additional five 
times, spanning Q a /(Q a +Q w ) values from 0.0 through nearly 0.5. The results of these 
computations are compiled in Figure 41 and Figure 42. 
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Figure 38. TECPLOT plot of residual history for two-phase TUTORIAL5145 simulation. 



Figure 39. TECPLOT plot of drag force history for two-phase TUTORIAL 5145 simulation. 
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iter fid ru rv rw rp 

ra 

rh rk 

re ... 

10000 1 5 . 809e-03 2.045e-03 2.639e-03 6.509e-01 

10000 2 3 . 713e-01 4.032e-01 1.596e-01 6.509e-01 

Awetted dot i: 2 . 236872 1 880023e-05 m A 2 

Viscous drag force: 2 . 7868256856116e+01 

Pressure drag force: 1 ,2050646663068e+01 

Total drag force: 3 . 9918903519184e+01 

1 . 4 82e-02 
1 . 099e-02 

0 . 000e+00 6 . 091e-04 
0.000e+00 0.000e+00 

2 . 577e-02 ... 

0 . 000e+00 

Figure 40. Snippets from standard output (n.out) for two phase TUTORIAL5145 simulation, illustrating case 

specific drag output. 



Figure 41. Surface pressure contours and a gas =0.5 isosurface for Uj„j=0.3 m/s TUTORIAL_5145 

simulation. 
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Figure 42. Predicted Total Hull Drag vs. Non-Dimensional Gas Injection for 5415 Simulations 



Figure 43. Predicted Total Drag Reduction vs. Non-Dimensional Gas Injection for 5415 

Simulations 
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Tutorial Case 4: Compressible Flow over a Bump 

This case demonstrates the application of NPHASE-PSU to a highly compressible 
flow field. Specifically, a calorically perfect air flow at an inlet Mach number of 0.7 flows 
inviscidly over a circular arc bump in a channel. Inviscid runs are accomplished by 
specifying Symmetry boundaries on walls and setting the viscosity to a tiny number. 

To start the tutorial the user needs to go to the TUTORIAL MACH BUMP directory. 
There are three files there when the software is unpacked: nphase.dat, cobalt. inp, and 
cobalt.be. The first step is to execute fump (as described in the Preprocessing section and 
Tutorial _1 section above), using 12 domains and a scaling factor of 1.0. This step generates 
grid+topology files, unphase. gridOOO - unphase. gridOl 1. 

To define initial conditions the following procedure was invoked: 

1) The keywords "use stagnation enthalpy as primitive energy 
transport variable" means just that, so stagnation enthalpy must be set at 
the inlet and as an initial guess for enthalpy. This is the standard approach used 
for highly compressible flow (but not necessarily for weakly compressible flow 
where static enthalpy can be more convenient [e.g., thermal driven convection]). 

2) Assuming an inlet static temperature and pressure of 293. 15K (20°C)and 
101325 Pa, we compute an inlet density from the perfect gas law and the 
defined specific heat ratio and gas constant for air (y= 1 .4, R=287 J/kg*K) 4 
p= 1.204328093 kg/m which appears in the inlet and pressure patch attributes. 

3) The inlet and initial velocity is computed taking an inlet Mach number of 0.7 
and computing the sound speed from the inlet static temperature for calorically 
perfect air (a=sqrt(yRT)). 

4) The stagnation enthalpy is computed by determining the static enthalpy for 
calorically perfect air from h=CpT, with Cp=yR/(y-l), and h0=h+u /2. 

The code is run on 12 processors for 2500 iterations by submitting an appropriately 
modified run.nphase script to the queue. The run converges fairly slowly due to the fine 
grid, 2 nd order convection numerics, the use of simple (i.e. non-characteristic) 
inflow/outflow boundary conditions, and the sharp shock that arises. 

The user can then run emerge (12 processor, 1 field) and migrate the merged data 
files (engold.geo, engold.uOO.Esca, engold.vOO.Esca, engold.wOO.Esca, 
engold.uvw00.Evec, engold.rOO.Esca, engold.pOO.Esca, engold.mach00.Esca, 
engold.h00.Esca, engold.tOO.Esca, and case file (engold.case) to the location he plans on 
postprocessing using ENSIGHT. 

Figure 44 shows a contour plot of the predicted Mach number contours. 
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#case title: 

#INVISCID BUMP 

iterations to perform 2500 

number of fields 1 

#initialize run with restart file 

produce ensight output 

restart file iteration write frequency 1000 

dont perform wall match logic 

#pamb=101325 Pa, Tamb=293.15 K 
inlet patch 1 0 

240.24145833 0. 0. 1.204328093 1. 0. 0. 0. 323327.154 0 0 #stagnation enthalpy 

pressure patch 1 0 

101325. 1.204328093 1. 0. 0. 323327.154 0. #stagnation enthalpy 
constant fluid molecular viscosity 1.0e-15 

use stagnation enthalpy as primitive energy transport variable 0 

strongly coupled compressibility 
perfect gas compressibility parameters 0 
287. 1.4 

single phase heating on 

continuity error treatment for enthalpy 
print umin and umax 
print tmin and tmax 

#clips on velocity -500. 500. -500. 500. -500. 500. 

spatialdiscretizationmomentum 2 
spatialdiscretizationenthalpy 2 

function entry/exit echo off 

#relaxation factor for p .3 
relaxation factor for u 0.7 
relaxation factor for v 0.7 

solver choice for velocity components jacobiuvw 
solver choice for pressure petsc 
#solver choice for pressure jacobi 

parallel strategy for pressure corrector: matrixlevel 

initialize u field 240.24145833 
initialize p field 101325. 

initialize h field 323327.154 #stagnation enthalpy 

Figure 44. nphase.dat file for TUT ORIAL M ACH BUMP. 
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Tutorial Case 5: Flow around a Shrouded Gear 


To demonstrate application of NPHASE-PSU to gear windage simulation, we have 
provided a tutorial case of a 72-tooth shrouded spur gear investigated by Diab et al. (2004). 
Considering the symmetry of the problem, only one side of a single tooth sector is 
simulated, decreasing simulation size substantially. The simulation is performed in a 
reference frame rotating with the gear, so that a stationary mesh can be used with a moving 
wall boundary condition applied on the shroud. Figure 46 shows the nphase.dat fde used 
for this case. Because the simulated flow is incompressible and there are no inlets or 
outlets, pressure must be set to a reference value (0) somewhere within the domain for the 
equations of motion to remain well-posed. We have chosen to anchor the pressure value 
implicitly (appearing in the pressure equation matrix). This is specified through inclusion 
of the line, anchor pressure implicitly. 

Attributes specific to gear windage simulation are specified in the line 
employ diab gear modeling 1000.0 2 72 o . The first argument is the gear rotation 
speed in radians per second. The second argument is an integer denoting the strategy used 
for the rotation. 0 is a stationary grid with a stationary reference frame, and implicit wall 
motion, but is applicable only to cylinders. 1 is for a stationary grid in the relative reference 
frame with absolute velocities. 2 is a stationary grid in the relative frame with relative 
velocities. 3 is a rotating grid in the absolute reference frame with absolute velocities. Here, 
we have chosen to use option 2, for a stationary grid in the relative frame with relative 
velocities. The third argument is the number of teeth. The fourth argument is a Boolean 
toggle denoting whether the output should be made relative to the first time-step. As this is 
a steady state simulation, this parameter is irrelevant and is set to zero. Due to the lack of 
inlets and outlets, the cyclic boundary conditions, and the manner in which flow is driven, 
this simulation takes many iterations to converge. To maintain numerical stability, it is best 
to run it with first order discretization for the first 50,000 iterations, then uncomment the 
second order discretization commands in nphase.dat, and run it again, re-initializing from 
restart files. 

To run this case, NPHASE requires that the high and low cyclic boundary faces be 
sorted. First, sort_cyclic is run, which generates cobalt. inp. sort. Next, this file must be re- 
named cobalt. inp. From here, fump is run. The current case is decomposed for 12 
processors using METIS and a scale of 1 . The sector angle for cyclic simulation is 5 
degrees for a 72-tooth gear, so this is input when prompted. After running fump to generate 
unphase. gridOOO - unphase. gridOl 1, the run script can be appropriately modified for the 
batch scheduling system and mpiexec used, and then submitted to the queue. After running 
with first order accuracy, nphase.dat should be modified for second order accuracy and 
initializing from restart files. To move the nphase_restart_out### files to 
nphase_restart_in###, run mrf. The run script can then be re-submitted to the queue to 
achieve a second order accurate solution. Figure 47 shows streamlines colored by pressure 
for this simulation. The streamlines were created in EnSight by applying computational 
symmetry to the solution as a post-processing step. 
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# Case title 

# Diab gear windage tutorial 

iterations to perform 50000 

restart file iteration write frequency 5000 

#initialize run with restart file 

number of fields 1 


produce ensight output 

employ diab gear modeling 1000.0 2 72 0 

#spatial discretization momentum 2 
#spatial discretization turbulence 2 


dont perform wall match logic 

anchor pressure implicitly 

pressure gradient at walls extrapolation 

constant fluid molecular viscosity l.e-5 
constant fluid density 1 . 

turbulence model for each field 1 

turbulence model reynolds number regime for each field 0 
production destruction ratio clip 500 


solver choice for velocity components jacobiuvw 
solver choice for pressure petsc 

parallel strategy for pressure corrector: matrixlevel 


#relaxation factor for u .7 
relaxation factor for k 0.7 
relaxation factor for e 0.7 


initialize 

initialize 

initialize 

initialize 


u field 0. 
p field 0. 
k field 0.015 
e field 5.0 


Figure 46. nphase.dat file for TUTORIALGEAR. 



P [Pa] 

1300 |=| 
175 
-950 
-2075 
-3200 L 


Figure 47. Streamlines colored by pressure from shrouded gear tutorial. 
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Other Items of Interest 


Running Two-Dimensional Problems 

NPHASE-PSU can be run in two-dimensional mode quite easily by specifying a 
one-element thick grid and defining symmetry boundaries on the front and back faces. It 
does not matter how thick the element is. So, for example, a pure 2D triangular mesh will 
extrude to a one-element thick prism mesh. 

Building User Specific/Case Specific Postprocessing 

Often an NPHASE-PSU user is interested in case specific output, or even case 
specific coding that modifies some element of execution (say grid motion, boundary 
condition, hard coded initialization, etc...). For these situations, this section documents the 
procedures to incorporate such coding in a fashion that is accessed from the front end, i.e., 
does not affect the execution of the code for cases where user supplied keywords are not 
supplied. The process is illustrated with an example - defining a user specific output for the 
5415 simulation carried out in the preceding tutorial. 

The first step is to define a new keyword and attendant integer pointer flag. In this 
example we add the following lines to the main NPHASE-PSU data C-structure in 
NPHASE-PSU, struct data, which is defined in nphase_struct.h: 

// 5415 specific stuff: 
int *imodel5415_modeling ; 

Next, storage is allocated for this pointer in constmemory.c: 

// 5415 specific: 

var.imodel5415_modeling= (int *) malloc(sizeof(int)) ; 

The value of this integer is initialized to 0 (i.e., not set) in set_parameters.c: 

int *imodel5415_modeling=var.imodel5415_modeling ; //type definition 


// 54 1 5 specific 
*imodel54 1 5_modeling=0; 

The value of this integer pointer is broadcast to all processors in broadcast, c: 

int *imodel5415_modeling=var.imodel5415_modeling ; //type definition 


MPI_Bcast(imodel54 1 5_modeling, 1 ,MPI_INT,0,MPI_COMM_WORLD); 

User access to setting this parameter is accomplished by defining a keyword in 
author * ***.c, where *** is either “abed”, “efgh”, “ijkl”, “mnop”, “qrstu”, or “vwxyz” 
depending on the fist letter of the keyword. Here we define the keyword “employ model 
5145 modeling” in author efgh.c: 
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int *imodel5415_modeling=var.imodel5415_modeling ; //type definition 


/* */ 

else if(stmcmp(keyword,"employmodel54 1 5modeling" 
,strlen(keyword)) = 0) { 
*imodel5415_modeling = 1; 

*idid = 1; 

fprintf(fo,"%s",line) ; 


Note that keyword is defined with no spaces or special characters. 

Now that all of the data structure and front end hooks are in place, the user can 
build or modify a function to execute what they wish. In this case, a new function is 
defined, model5415_output.c, the source of which appears in Figure 43. This function is 
called, in this case, at the end of every iteration, from nphase.c: 

int *imodel5415_modeling=var.imodel5415_modeling ; //type definition 


if(*imodel54 1 5_modeling= 1 )model54 1 5_output(); 

Once these coding modifications have been made, the object model5415_output.o is 
add to the file Obj nphase, and the entire code is recompiled using make -f makefile 
nphase (options). The 5415 tutorial presented above includes the use of this coding as 
shown in the nphase.dat files included in Figure 34. The output generated by this coding is 
shown in Figure 36 and Figure 39. 
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#include "nphase_struct.h" 

#include "mpi.h" 

#include <math.h> 

#include <stdio.h> 

extern struct boundary_patch wall; 
extern struct boundary_patch porwall ; 
extern struct data var; 
extern int myid; 

int model5415_output ( ) 

{ 

/* 

called from; 
nphase 

V2 . 0 baseline code 
initials comment 

rfk 

*/ 

int *nf ield=var . nf ield; 
int *nnode=var . nnode ; 
int ibf, nstride, inode, fstride, if ield; 
real sumforce, sumforceO, sumflux, sumfluxO; 
real awetted, awettedO ; 
real rinf , uinf , cd; 
real *a=var.a; 
real *u=var.u; 

entered ( "model5415_output" ) ; 

uinf=2 . 2134 ; 
rinf=1000; 

/ / wetted area 

awetted=0 . ; 

for (ibf=0; ibf<=wall . nbcf ace-1 ; ++ibf) { 
inode=wall . bcf ace_n [ ibf] ; 
awetted+=wall . bamag [ ibf ] ; 

} 

for ( ibf =0 ; ibf<=porwall . nbcf ace- 1 ; ++ibf ) { 
inode=porwall .bcface_n [ibf] ; 
awetted+=porwall . bamag [ ibf ] ; 

} 

MPI_Reduce (&awetted, SawettedO, 1 , mpireal, MPI_SUM, 0,MPI_COMM_WORLD) ; 
if (myid==0) printf ( "Awetted: %20.13e m A 2\n" , awettedO) ; 

// drag force on boat 

sumforce=0 . ; 

for ( if ield=0 ; if ield<=*nf ield-1 ; ++i field) { 
fstride=if ield*wall . nbcf ace ; 
nstride=if ield* *nnode; 

for ( ibf =0 ; ibf<=wall . nbcface-1 ; ++ibf ) { 
inode=wall . bcf ace_n [ ibf ] ; 

sumforce+=wall . tmlt [ ibf +f stride] * * (a+inode+nstride) * wall . bamag [ ibf ] * * (u+inode+nstride) ; 

} 

for (ibf=0; ibf<=porwall . nbcf ace-1 ; ++ibf) { 
inode=porwall .bcface_n [ibf ] ; 

sumforce+=porwall . tmlt [ibf+f stride] * * (a+inode+nstride) * porwall .bamag [ibf ] * * (u+inode+nstride) ; 

} 

} 

MPI_Reduce (&sumforce, &sumforce0, 1 , mpireal, MPI_SUM, 0 , MPI_COMM_WORLD) ; 
if (myid==0 ) printf ( "Drag force: %20 . 13e\n" , sumforceO ) ; 
if (myid==0) cd=sumforceO/ ( . 5*rinf *uinf *uinf *awetted0 ) ; 
if (myid==0) printf ("CD=: %20 . 13e\n", cd) ; 

exiting ( "model5415_output" ) ; 

return 0 ; 

} 

Figure 48. model5415_output.c 
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T urbomachinery 

Documentation not yet written. 


Homogeneous Multiphase Flow 

Documentation not yet written. 
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NPHASE-PSU Reference Manual 


Control Commands 

This section includes a description of most of the keyword commands available in 
NPHASE. These appear in nphase.dat in free format. They need not appear in any 
particular sequence. It is not exhaustive and many commands are not likely to be used by 
most users. 

1.1. Comment Card 

Any line in the “nphase.dat” input file that starts with a “#” is considered a comment card. The remainder of 
the line is ignored. 

1.2. Job Control 

1.2.1. Case Title 

Specifies a title for the output of the job. Optional. A maximum of 132 characters can be used in the title. The 
title itself must appear on following line 

case title: 

$ Title 


1.2.2. Iterations to Perform 

Specifies the number of iterations to perform in job. For a restart job, the specified number of iterations will 
be performed after the restart file is read. For a time accurate simulation, this is the number of inner iterations 
per physical timestep. 

iterations to perfonn $NITER 

1.2.3. Over Write Boundary Conditions on Restart 

During a restart from a previous run (See “initialize run with restart file” command) the boundary conditions 

are specified by the information in the restart file. If the user wants the boundary conditions specified by the 

“nphase.dat” input file, the over write boundary condition on restart command must be used. The command 

can be used on various boundary conditions as specified below. 

over write inlet patch boundary conditions on restart 

over write far field patch boundary conditions on restart 

over write pressure patch boundary conditions on restart 

over write wall patch boundary conditions on restart 

over write porous patch boundary conditions on restart 

1.2.4. Number of Physical Time Steps 

The number of time steps for the case can be specified with this command. 

number of physical time steps $NTSTEPS 

1.2.5. Physical Time Step in Seconds 

This command is used to specify a constant time step size for a transient analysis. 

physical time step in seconds $DT 

1.2.6. Time Accurate Simulation 

Specifies that a transient analysis is to be run. The number of time steps and time step size can be specified 
with the “number of physical time steps” and “physical time step size in seconds” command 
time accurate simulation 
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1.2.7. Transient File Write Frequency 


Specifies the number of time steps between saving a time step numbered restart file. This command can be 
used to save restart files for postprocessing animation. 


restart file write frequency $NRST FREQ 

1.2.8. Volume Fraction Normalization 

Specifies normalization of volume fraction equation. Normalization will ensure the sum of the volume 
fraction of all fields equals one. 

carver volume fraction normalization 

Specifies that volume fraction normalization is not used. 

do not employ carver volume fraction normalization 

1.2.9. Continuity Error Treatment for Momentum 

Adds terms to the LHS and RHS of the discrete momentum equations that ensures that if the linear solver is 
brought in then no spurious sources of momentum will arise due to mass imbalances 

continuity error treatment for momentum 

1.2.10. Continuity Error Treatment for Enthalpy 

Adds terms to the LHS and RHS of the discrete enthalpy equations that ensures that if the linear solver is 
brought in then no spurious sources of enthalpy will arise due to mass imbalances 

continuity error treatment for enthalpy 

1.2.11. Momentum Cross Diffusion 

This command specifies that the nonorthogonal components of the viscous shear term are included in the 
momentum equation. For a general grid, this command should be included in the input deck. See equation 
Error! Reference source not found, in the theory manual. 

cross diffusion included in momentum equations 

1.2.12. Number of PISO Correction Steps 

The number of PISO corrections steps is required when the PISO algroithm is employed (see piso employed 
command). The PISO algorithm is a more recent varient of the SIMPLE algroithm. 

number of piso corrections $NPISO STEPS 

1.2.13. Parallel Strategy for Pressure Correction 

Specifies whether the pressure correction equation is solved implicitly across all inter-processor boundaries. 

parallel strategy for pressure correction: matrix level 


parallel strategy for pressure correction: explicit partition boundary update 

1.2.14. PISO Algorithm Employed 

The PISO algorithm is a more recent variant of the SIMPLE algorithm. The PISO algorithm can be specified 
using this command. The number of PISO correction steps can be specified using the “number of piso 
corrections” command. 

piso employed 


1.2.15. Relaxation Factor 

The relaxation factor command is used to improve numerical convergence by the addition of numerical 
damping to eliminate oscillations and improve stability in the solution. The form of the relaxation factor 
command is: 

relaxation factor for u $RFU 
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relaxation factor for v $RFV 

relaxation factor for w $RFW 

relaxation factor for a $RFA 

relaxation factor for k $RFK 

relaxation factor for e $RFE 

relaxation factor for h $RFH 

Values for the false time step command are specified as: 

$RFU - User input value for false time step in “u” momentum equation 
$RFV - User input value for false time step in “v” momentum equation 
$RFW - User input value for false time step in “w” momentum equation 
$RFA - User input value for false time step in phasic volume fraction equation 
$RFK - User input value for false time step in turbulent kinetic energy equation 
$RFE - User input value for false time step in turbulent dissipation rate equation 

1.2.16. SIMPLEC Factor 

The use of this command can modify the numerical method to use the original SIMPLE algorithm 
($SIMPLEC=0.) or a varient call the SIMPLEC algorithm ($SIMPLEC=1.). The SIMPLEC algorithm is the 
default value. 


simplec factor $ SIMPLEC 

1.2.17. Solver Choice for All Scalars Jacobi 

This command specifies the use of the Jacobi algorithm for solving the volume fraction, k-e turbulence and 
energy equations. This is the default solver for the volume fraction, k-e turbulence and energy equations. 

solver choice for all scalars jacobi 

1.2.18. Solver Choice for Enthalpy Jacobi 

This command specifies the use of the Jacobi algorithm for solving the energy equation. This is the default 
solver for the energy equation. 

solver choice for all scalars jacobi 

1.2.19. Solver Choice for Pressure Jacobi 

This command specifies the use of the Jacobi algorithm for solving the pressure correction equations in the 
SIMPLE algorithm. This is the default solver for the pressure correction equations 

solver choice for pressure jacobi 

1.2.20. Solver Choice for Turbulence Scalars Jacobi 

This command specifies the use of the Jacobi algorithm for solving the k-e turbulence equations. This is the 
default solver for the k-e turbulence equations 

solver choice for turbulence scalars jacobi 

1.2.21. Solver Choice for Velocity Components Jacobi 

This command specifies the use of the Jacobi algorithm for solving the momentum equations in the SIMPLE 
algorithm. This is the default solver for the momentum equations 

solver choice for velocity components jacobi 

1.2.22. Solver Sweeps 

This command is used to specify the number of linear solver sweeps for each linear equation solver. 
Command are available to specify the number of solver sweeps separately for each equation or in 
combination.. 

solver sweeps for u $NSWEEP field 1, $NS WEEP field 1,...,$NS WEEP fieldN 
solver sweeps for v $NSWEEP field 1, $NS WEEP field 1,...,$NS WEEP fieldN 
solver sweeps for w $NSWEEP field 1, $NS WEEP field 1,...,$NS WEEP fieldN 


NASA/CR— 2014-2 1665 1 


77 

119 


solver sweeps for p $NSWEEP_fieldl, $NSWEEP_fieldl,...,$NSWEEP_fieldN 

solver sweeps for a $NSWEEP_fieldl, $NSWEEP_fieldl,...,$NSWEEP_fieldN 

solver sweeps for h $NSWEEP_fieldl, $NSWEEP_fieldl,...,$NSWEEP_fieldN 

solver sweeps for tke $NSWEEP_fieldl, $NSWEEP_fieldl,...,$NSWEEP_fieldN 

solver sweeps for tds $NSWEEP field 1, $NS WEEP field 1,...,$N SWEEP fieldN 

1.2.23. Spatial Discretization 

This command is used to specify the spatial discretization used for the convective term in the governing 
equations. The default discretization is the hybrid scheme. 

spatial discretization momentum $ISPATIAL 
spatial discretization turbulence $ISPATIAL 
spatial discretization volume fraction $ISPATIAL 

spatial discretization turbulence $ISPATIAL 

Values for the spatial discretization command are specified as: 

$ISPATIAL = 0 ( Use 1st order hybrid scheme) 
l(Use 1st order upwind scheme) 

2(Use 2nd order upwind scheme) 

1.2.24. Temporal Discretization 

This command is used to specify the spatial discretization used to transfrom the convective term in the 
governing equations. The default discretization is the hybrid scheme. 

temporal discretization momentum $ITEMPOR 
temporal discretization turbulence $ITEMPOR 
temporal discretization volume fraction $ITEMPOR 
temporal discretization turbulence $ITEMPOR 

Values for the temporal discretization command are specified as: 

$ITEMPOR = 1 (Use 1 st order [Euler Implicit] scheme - default) 

$ITEMPOR = 2 (Use 2 nd order backward difference in time) 

1.2.25. Thin Layer Approximation Employed/ Not Employed 

Specifies whether or not a thin layer approximation is employed in the construction of viscous terms, 
thin layer approximation employed 
thin layer approximation not employed 


1.3. Geometry Control 

1.3.1. Cylindrical Coordinates 

Specifies use of cylindrical coordinates (R-Z). Where axial flow (Z) is in the (X) coordinate direction, and 
radial (R) flow is in the (Y) coordinate direction. Requires (Y=0) to align with (R=0). 

Cylindrical Coordinates 


1.4. Boundary Conditions 

1.4.1. Cyclic Boundary 

All cyclic boundaries are specified in the “unphase. grid” file. 

1.4.2. Far Field Boundary 

1.4.3. Far Field Patch 

All far field boundaries are specified in the “unphase. grid” file. 

farfield patch $NCARD, $FACEID 
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$U, $V, $W, $RHO, $VF, $TKE, $TDS, $H,$P // for field one, repeat NCARD times for other 
fields 

Specified far field boundary conditions for all farfield boundary faces with FACEID. Far field values are 
specified as: 

$NCARD - Number of input lines imediately following. Usually equals NFIELD. 

SFACEID - Identification of particular far field patch, (e.g., FACEID = 3 for Farfield 03) 

$U - Cartesian “u” velocity at far field boundary (m/s) 

$V - Cartesian “v” velocity at far field boundary (m/s) 

$W - Cartesian “w” velocity at far field boundary (m/s) 

$RHO - Fluid Density at inlet (kg/m 3 ), used for backflow conditions only 

$VF - volume fraction at far field boundary 

$TKE - Turbulent Kinetic Energy at far field boundary (m 2 /s 2 ) 

STDS - Turbulence Dissipation at far field boundary 
$H - Enthalpy at far field boundary (deg-K) 

$P - Pressure at far field boundary (Pa) 

1.4.4. Inlet Boundary 
1.4.4.1.Inlet Patch 

The Inlet Patch command is used to specify uniform inlet boundary conditions. All inlet boundaries are 
specified in the “unphase. grid” file. 

Inlet Patch $NCARD, $FACEID 

$U, $V, $W, $RHO, $VF, $TKE, $TDS, $H, $P, $T1, $T2 // for field one, repeat NCARD 
times for other fields 

Specified inlet flow boundary conditions for all inlet boundary faces with FACEID. Inflow values are 
specified as: 

SNCARD - Number of input lines imediately following. Usually equals NFIELD. 

SFACEID - Identification of Faces to apply inlet patch. Face ID is specifed in “unphase. grid” file. 

$U - Cartesian “u” velocity at inlet (m/s) 

$V - Cartesian “v” velocity at inlet (m/s) 

$W - Cartesian “w” velocity at inlet (m/s) 

$RHO - Fluid Density at inlet (kg/m 3 ), used for backflow conditions only 

$VF - volume fraction at inlet 

$TKE - Turbulent Kinetic Energy at inlet (m 2 /s 2 ) 

STDS - Turbulence Dissipation at inlet 

$H - Enthalpy at inlet (deg-K). Specified either stagnation or static depending on ienergy variable 

$P - Stagnation pressure at inlet 

$T1 =0,1 for subsonic, supersonic treatment 

$T2 = 0,1 for incompressible, compressible 

1.4.4.2.Inlet Velocities Specified in Absolute Frame. 

inlet velocities specified in absolute frame 

1.4.4.3.Inlet Velocity Specified in Cylindrical Coordinates. 

inlet velocities specified in cylindrical coordinates 

1.4.4.4.Inlet Velocities Specified Transiently. 

inlet velocities specified transiently 

1.4.4.5.Inlet Profile Specified in inflow.pro 

inlet profile specified in inflow.pro 

1.4.4.6.Inlet Profiles 2D specified in inflow.pro 

inlet profile 2d specified in inflow.pro 
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1.4.5. Porous Boundary 

1.4.5.1. Porous Wall Patch 

The Porous Wall Patch command is used to specify uniform boundary conditions through a porous surface. 

All porous boundaries are specified in the “unphase. grid” file. 

porous wall patch $NCARD, $FACEID 

$U, $V, $W, $RHO, $VF, $TKE, $TDS, $T, $PERMA // for field one, repeat NCARD times 
for other fields 

Specified inlet flow boundary conditions for all inlet boundary faces with FACE1D. Inflow values are 
specified as: 

$NCARD - Number of input lines imediately following. Usually equals NF1ELD. 

SFACEID - Identification of Faces to apply inlet patch. Face ID is specifed in “unphase. grid” file. 

$U - Cartesian “u” velocity through porous wall (m/s) 

$V - Cartesian “v” velocity through porous wall(m/s) 

$W - Cartesian “w” velocity through porous wall (m/s) 

$RHO - Fluid Density through porous wall (kg/m 3 ) 

$VF - volume fraction at through porous wall 

$TKE - Turbulent Kinetic Energy through porous wall(m 2 /s 2 ) 

STDS - Turbulence Dissipation through porous wall 
$T - Temperature through porous wall (deg-K) 

SPERMA - permeability of porous wall 

1.4.6. Pressure Boundary 

1.4.6.1. Pressure Patch 

The Pressure Patch command is used to specify uniform pressure boundary conditions. All pressure 
boundaries are specified in the “unphase. grid” file. 

Pressure Patch $NCARD, $FACEID 

$P, $RHO, $VF, $TKE, $TDS, $H, $P0, $TYPE // for field one, repeat NCARD times for 
other fields 

Specified inlet flow boundary conditions for a patch of faces. Inflow values are specified as: 

$NCARD - Number of input lines imediately following. Usually equals NF1ELD. 

SFACEID - Identification of Faces to apply inlet patch. Face ID is specifed in “unphase, grid” file. 

$P - Pressure at pressure boundary (Pa) 

$RHO - Fluid Density at pressure boundary (kg/m 3 ), used for backflow conditions only 

$VF - volume fraction at pressure boundary, used for backflow conditions only 

$TKE - Turbulent Kinetic Energy at pressure boundary (m 2 /s 2 ), used for backflow conditions only 

STDS - Turbulence Dissipation at pressure boundary, used for backflow conditions only 

$H - Enthalpy at pressure boundary (J/kgK), used for backflow conditions only 

$P0 - Stagnation pressure at pressure boundary (Pa), currently unused 

$TYPE - Included for future support. Just use 0 for now 

1.4.7. Symmetry Boundary 

All symmetry boundaries are specified in the “unphase, grid” file. 

1.4.8. Wall Boundary 

1.4.8.1. Wall Patch 

The Wall Patch command is used to modify wall boundary conditions. All wall boundaries are specified in 
the “unphase. grid” file. 

Wall Patch $NCARD, $FACEID 

$XTW, $XMW, $XFW, $XFW0, $XFW1, $XFW2 //heat transfer card 
$XTW_STRUCT, $XFW_TENSION //wall structural parameters card 
$XTW VEL, $XFW XVEL, $XFW VVEL, $XFW ZVEL //wall motion card 
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//for field one, repeat these 3 lines NCARD times 

for other fields 

NOTE: one does not need to specify any wall attributes if the wall is adiabatic, rigid and statinoary. If any 

of these attributes is non-default, all must be specified 

Specified wall boundary conditions for all wall boundary faces with FACEID. Wall attributes are specified 
as: 

$NCARD - Number of input lines imediately following. Usually equals NFIELD. 

SFACEID - Identification of Faces to apply wall patch. Face ID is specifed in “unphase, grid” file. 

$XTW - Type of wall boundary: 

iwt=0: adiabatic 

iwt=l: specify temperature 

iwt=2: specify heat flux (xfw= 0. = adiabatic) 

iwt=3 : specify heat transfer coeficient (not implemented yet) 

iwt=4: conjugate heat transfer across thin matching walls (continuous flux on either side of wall) 

iwt=5: conjugate ID heat transfer across thickness of assumed material adjacent to wall face with temp 

specified on outside face of that material 

iwt=6: conjugate ID heat transfer across thickness of assumed material adjacent to wall face with htc and 
tfilm specified on outside face of that material 

$XMW - If this is a matching wall this integer is the the negative of the faceid of the matching wall. If this is 
not a matching wall it is ignored. 

$XFW - Real # parameter: 

if iwt=l, xfw=specified temperature in °K 

if iwt=2, xfw=specified heat flux in W/m 2 

if iwt=3, xfw=specified heat transfer coefficient (not implemented) 

if iwt=4, xfw is unused since heat flux is continuous across matching internal walls 

if iwt=5, xfw=specified temperature on outside of material assumed adjacent to wall patch 

if iwt=6, xfw=specified film temperature (K) on outside of material assumed adjacent to wall patch 

$XFW0 - unused for iwt=0->4 

if iwt=5, xfwO=shell thermal conductivity (J/m*s*K) of material assumed adjacent to wall patch 
if iwt=6, xfw0= heat transfer coefficient (W/m 2 *K) on outside of material assumed adjacent to wall patch 
$XFW1 - unused for iwt=0->4 

if iwt=5, xfwl=shell thermal conductivity (J/m*s*K) of material assumed adjacent to wall patch 
if iwt=6, xfwl= shell thickness (m) of material assumed adjacent to wall patch 
$XFW2 - unused for iwt=0->5 

ifiwt=6, xfw2=shell thermal conductivity (J/m*s*K) of material assumed adjacent to wall patch 

1.5. Initial Conditions 

1.5.1. Initialize with Restart 

The restart file “nphase restart in” will be used to specify the boundary and interior node values, 

initialize run with restart file 

1.5.2. Initial u Velocity 

This command initializes the u velocity to a constant value specified by the user:. 

initialize u velocity $u_fieldO, $u_fieldl,...,$u_fieldN 

1.5.3. Initial v Velocity 

This command initializes the v velocity to a constant value specified by the user:. 

initialize u velocity $v_field0, $v_fieldl,...,$v_fieldN 

1.5.4. Initial w Velocity 

This command initializes the w velocity to a constant value specified by the user:. 

initialize u velocity $w_fieldO, $w_fieldl,...,$w_fieldN 
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1.5.5. Initial Pressure Field 

This command initializes the pressure field to a constant value specified by the user:. 


initialize p field $w fieldO, $w fieldl,...,$w fieldN 

1.5.6. Initial Volume Fraction Field 

This command initializes the volume fraction field to a constant value specified by the user:. 

initialize volume fraction field $VF fieldO, $VF fieldl,...,$VF fieldN 

1.5.7. Initial Interfacial Area Density 

This command initializes the volume fraction field to a constant value specified by the user:. 

initialize ai field $AI fieldO, $AI fieldl,...,$AI fieldN 

1.5.8. Initial Turbulent Kinetic Energy (tke) Field 

This command initializes the tke field to a constant value specified by the user:. 

initialize tke field $TKE fieldO, $ TKE fieldl,...,$TKE fieldN 

1.5.9. Initial Turbulent Dissipation Rate (tds) Field 

This command initializes the turbulent dissipation rate field to a constant value specified by the user:. 

initialize tds field $TDS fieldO, $TDS fieldl,...,$TDS fieldN 

1.5.10. Initial Enthalpy Field 

This command initializes the enthalpy field to a constant value specified by the user:. 

initialize h field $H fieldO, $H fieldl,...,$H fieldN 

1.5.11. Hard Coded Initialization 

hard coded initialization 


1.6. Physical Models 

1.6.1. Environmental Properties 

1. 6.1.1. Laminar Flow 

This command specifies that the flow is to be treated as laminar flow., 

laminar flow 

1.6.1.2. Gravity Vector 

This command specify the direction and magnitude of gravity. 

gravity vector $gx, $gy, $gz 

Where, 

$gx == “x” direction gravity magnitude, 

$gy == “y” direction gravity magnitude, 

$gz == “z” direction gravity magnitude 

1. 6.1.3. Employ Modified Pressure 

This command specifies that a conventional gravity modified pressure treatment will be used. 

employ modified pressure 

1.6.1.4. Reference Density for Modified Pressure 

This command specifies the reference density for the modified pressure treatment. 

reference density for modified pressure $rhoref 
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I.6.I.5. Employ Simple Hydrostatic Pressure Treatment 


This command specifies that a simple single phase hydrostatic initialization of the pressure will be used based 
on local coordinate (Vp=pgx) and modified pressure reference density. 


employ simple hydrostatic presure treatment 

1.6.1. 6. System Rotation About X Axis in Radians Per Second 

Specifies the system rotation about the x-axis in radians per second. 

system rotation about x axias in radians per second $XRADS 

1.6.2. Fluid and Solid Properties 

1.6.2. l.Constant Fliud Molecular Viscosity 

The constant fluid molecular viscosity is used to specify the molecular viscosity for each field. 

constant fluid molecular viscosity $MU field 1, $ MU fieldl,...,$MU fieldN 

1.6.2.2.Constant Fluid Density 

The constant fluid density command is used to specify the fluid density for each field.. 

constant fluid density $RHO field 1, $RHO fieldl,.. 

1.6.2.3.Constant Fluid Surface Tension 

The constant fluid surface tension command is used to specify the fluid surface tension for each field.. 

constant fluid surface tension $SIGMA fieldl, $ SIGMA fieldl,. ..,$SIGMA fieldN 

1.6.2.4.Constant Fluid Specific Heat at Constant Pressure 

The constant fluid specific heat at constant pressure command is used to specify the fluid specific heat (Cp) 
for each field.. 

constant fluid specific heat at constant pressure $CP fieldl, $CP fieldl, ...,$CP fieldN 

1.6.2.5.Isothermal Compressibility Parameters 

This command is used to specify compressibility parameters. 

isothermal compressibility parameters $NCARD 

l // for field one, repeat NCARD times for other fields 

1.6.2.6.Solid Region Density 

The density of solid regions can be specified by this command. 

solid region density $RHO solids 

1.6.2.7.Solid Region Thermal Conductivity 

This command specifies the thermal conductive for solid regions. 

solid region thermal conductivity $K solids 

1.6.2.8.Solid Region Specific Heat 

This command specifies the specific heat for solid region. 

solid region specific heat $CP solids 

1.6.3. Turbulence Models 

1.6.3. l.Enforce Production Equals Dissipation 

enforce production equals dissipation 

1.6.3.2.Turbulent Flow High Reynolds Number k-s Turbulence model 

This command specified the use of the industry standard k-s turbulence model. 

turbulent flow high reynolds number k epsilon 
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1.6.3.3.Turbulent Flow Low Reynolds Number k-£ Turbulence model 


This command specified the use of a low Reynolds number version of the k-s turbulence model (Chien). 


turbulent flow low reynolds number k epsilon 

1.6.3.4.TurbuIent Flow High Reynolds Number q-co Turbulence model 

This command specified the use of the q- to turbulence model. 

turbulent flow high reynolds number q omega 

1.6.3.5.Turbulent Flow Low Reynolds Number q- a Turbulence model 

This command specified the use of a low Reynolds number version of the q-oi turbulence model (Coakley). 

turbulent flow low reynolds number q omega 

1.6.4. Multiphase Models 

I.6.4.I. Constant Field Characteristic Diameter 

The constant field characteristic diameter command is used to specify a single effective diameter for each 
field. This diameter is used in multiphase models such as drag, wall force and other models. It is ignored 
except for initialization when coalescence or breakup are implemented or if interfacial area transport is 
employed 

constant field characteristic diameter $D field 1, $D field 1,...,$D fieldN 

NOTE: some fields (i.e., continuous liquid field) may not use the characteristic diameter, but an input is 
required. 

1.6.4.2.Number of Fields 

This command is used to specify the total number of fields. In a two-fluid ensembled averaged model, the 
fields can either represent different phases, of various forms of the same phase (i.e., small bubbles, large 
bubbles, taylor bubbles or continuous vapor). 

number of fields $NFIELDS 

1.6.4.3.Interfacial Drag Force 

The interfacial drag force command is used to specify either a user defined or built-in model between any 
number of fields. At least one interfacial drag model is typically used for each field, however, none are 
requird and more than one is allowed. The form of the interfacial drag command is: 

interfacial drag model $NCARD 

$FIELDA,$FIELDB,$MODELID,$USERMULTIPLIER // for first drag model, repeat NCARD times 
for other drag models 


The input for the drag force command is defined as: 

$NCARDS - Number of input lines immediately following command line. Each line will specify a drag 
relation between fields. 

$FIELDA - First Field for drag model 
$FIELDB -Second field for drag model 
SMODELID - Specifies Drag Model 

0 - User Defined Drag Model (Use “user drag.c” routine to define drag modeling) 

1 - standard bubbly flow drag, C D =0.44 

2 - standard bubbly flow drag, C D =0.44 with virtual mass force intergrated 

3 - particle drag (solid spehere) 

4 - seawater bubble drag 

5 - fresh water bubble drag 

6 - contaminated fresh water bubble drag 

$USERMULTIPL1ER - User specified multiplier (To modify drag coefficient, default=1.0) 

I.6.4.4. Interfacial Non Drag Force 

The interfacial nondrag force command is used to specify either a user defined or built-in model between any 
number of fields. The form of the interfacial nondrag command is: 
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interfacial nondrag model $NCARD 

$FIELDA,$FIELDB,$MODELID,$MODELSUBID,$USERMULTIPLIER // for first nondrag 

model, repeat NCARD times for other nondrag models 
The input for the nondrag force command is defined as: 

$NCARDS - Number of input lines immediately following command line. Each line will specify a drag 
relation between fields. 

$F1ELDA - First Field for drag model 
$F1ELDB -Second field for drag model 
$MODELID - Specifies Non Drag Force Model 

0 - user specified model 

1 - Lift Force Non-drag Model 

2 - Volume Fraction Dispersion Non-drag Model 

3 - Wall Force Non-drag Model 

4 - Turbulence Non-drag Model 

5 - Virtual Mass Force Non-drag Model 
SMODELSUBID - Specifies subid for Non Drag Force Model 

SUSERMULTIPLIER - User specified multiplier (To modify non-drag coefficient, default=1.0) 

1.6.4.5. Wall Shear Apportionment Model 

The wall shear apportionment command is used to specify the fraction of wall shear associated with each 
field. Since, in many multiphase applications, the wall adjacent nodes may contain multiple fields, the fields 
in contact with the wall should be assigned the wall shear. The default for the code is the wall shear is 
apportioned by the local volume fraction in the wall adjacent nodes. In many applications (i.e., bubble flows, 
annular flows, ...) the user may assign all or a fraction of the wall shear to a particular field (or fields) with 
this command. 

wall shear apportionment model $FRACT_fieldl, $FRACT_field2,...,$FRACT_fieldN 
Where, $FRACT_field is the fraction of wall shear associated with each field. The sum of all $FRACT_field 
must add to one inorder to account for the wall shear. 

1.6.4.6. Bubble Cluster Drag Modification 

Specifies that the drag coefficient is modified according to equation Error! Reference source not found.. 

Bubble cluster drag modification 

1.6.4.7.Interfacial Area Transport Model for Each Field 

Specifies the interfacial area transport model for each field. 

interfacial area transport model for each field $iaint_fieldO, $iaint_fieldl, ..., $iaint_fieldN 

1.6.4.8. Interfacial Area Coalescence Model for Each Field 

Specifies the interfacial area coalescence model for each field. 

interfacial area coalescence model for each field $coalescence_aint_model_fieldO, 
$coalescence_aint_model_fieldl , . . . , $coalescence_aint_model_fieldN 

1.6.4.9. Interfacial Area Breakup Model for Each Field 

Specifies the interfacial area breakup model for each field. 

interfacial area breakup model for each field $breakup_aint_model_fieldO, 
$breakup_aint_model_fieldl , . . . , $breakup_aint_model_fieldN 


1.6.4.10. Initialize Interfacial Area Using Volume Fraction and 
Characteristic Length 

Specifies that the initial interfacial area density is defined using the local field volume fraction and 

characteristic length, not by “initialize ai field” 

initialize interfacial area using volume fraction and characteristic length 
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1.6.4.11. Interfacial Area Feeds Back Into Bubble Diameter Frequency 

Specifies the iteration frequency with which field bubble diameters are recomputed based on interfacial area 

density. 

interfacial area feeds back into bubble diameter frequency $aint_db_feedback 

1.6.4.12. Employ Rhie Chow for Dispersion Terms 

Specifies that the Rhie-Chow artificial dissipation parameter Error! Reference source not found, includes 
the turbulence dispersion model. 

employ rhie chow for dispersion terms 

1.6.4.13. 


Heat Transfer Models 

1.6.4.14. Single Phase Heating 

The use of this command will include the solution of the energy equation in the run. This command will 
allow only single phase heating (i.e., no mass transfer) in the analysis. 

single phase heating 

1.6.4.15. Adiabatic Flow 

The adiabatic flow command is used to specify mass and momentum conservation are only needed for the 

NPHASE code. 

adiabatic flow 

1.6.4.16. Constant Fluid Reference Enthalpy 

The constant fluid reference enthalpy command is used to specify a reference enthalpy for use in heat transfer 
models.. 

constant fluid reference enthalpy $HREF_fieldl, $HREF_freldl,...,$HREF_fieldN 

1.6.4.17. Constant Fluid Reference Temperature 

The constant fluid reference temperature command is used to specify a reference temperature for use in heat 
transfer models. 

constant fluid reference temperature $TREF_fieldl, $TREF_fieldl,...,$TREF_fieldN 

1.6.5. Mass Transfer Models 

1.6.5.1.Breakup sink for carrier field turbulence 

breakup sink for carrier field turbulence 

1.6.5.2.InterField Mass Transfer Models 

This command can be used to specify a mass transfer model. The form of the interfacial drag command is: 

interfield mass transfer models $NCARD 

$FIELDA,$FIELDB,$CARRIER,$MODELID,$USERMULTIPLIER // for first model, repeat 

NCARD times for other models 

The input for the drag force command is defined as: 

$NCARDS - Number of input lines immediately following command line. Each line will specify a drag 
relation between fields. 

$F1ELDA - First Field for mass transfer model 
$F1ELDB - Second field for mass transfer model 
SCARR1ER - Carrier field for mass transfer 
$M0DEL1D - Specifies Drag Model 
1 - Standard Mass Transfer Model 

SUSERMULTIPLIER - User specified multiplier (To modify mass transfer coefficient, default=1.0) 
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1.6.5.3.Mass Diffusion Coefficient Laminar 

mass diffusion coefficient laminar $MCOEF_fieldl, $MCOEF_fieldl,...,$MCOEF_fieldN 

1.6.5.4.Mass Diffusion coefficient Turbulent. 

mass diffusion coefficient turbulent $MCOEF_fieldl, $MCOEF_fieldl,...,$MCOEF_fieldN 

1.7. Additional Commands 


allow deforming walls 


anchor pressure 


ausm artificial dissipation 


boussinesq heating beta rhoOtO $BETA BOUSS, $RHO REF BOUSS, $T REF BOUSS 


build patran neutral file 


clip species fraction stolie between zero and one inclusive 


do not clip species fraction stolie between zero and one inclusive 


clip volume fraction stolie between zero and one inclusive 


do no tclip volume fraction stolie between zero and one inclusive 


compute wall proximities 


dont compute wall proximities 


bubbly drag model $IP1, $IP2 


constant fluid shear modulus $SFIRMOD_fieldl, $SHRMOD_fieldl,...,$SFIRMOD_fieldN 


constant fluid reference density $RHOREF_fieldl, $RHOREF_fieldl,...,$RHOREF_fieldN 


dimensionality offield view output $FV_OUT_DIM 


continuous fieldO turbulence drives all other fields 


des modifications 
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dont perform wall match logic 


drag and mass transfer in pseudo time step specification 


drag and mass transfer not in pseudo time step specification 


do not employ rhie chow for dispersion terms 


continuity error treatment for turbulence 


continuity error treatment for volume fraction 


continuity error treatment for species 


deforming region specification mid density modulus poisson $IMID, 

$L REG DEFORM, $R REG DEFORM, $E REG DEFORM, $N REG DEFORM 


deforming wall specification face id mid modulus poisson thickness $IMID, 
$L_WALL_DEFORM, $R_WALL_DEFORM, $E_WALL_DEFORM, 

$N WALL DEFORM 


defonning grid iteration frequency $ITER_FREQ 


des model control fsst flag $FSST_FLAG 


des model control c des $CDES 


des model control t scale des $TSCALE DES 


des model control sigma x des $S1GMAX_DES 


des model control chi des $CH1 DES 


des model control ch2 des $CH2 DES 


des model control ch3 des $CH3 DES 


breakup model mu $BREAKUP_MULT 


coalescence model mu $C MULT 
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coalescence model surfactant $C SURFACTANT 


coalescence model initial film thi$TFILM 0 


coalescence model final film thi$TFILM FINAL 


cfl number mear$CFL 


cfl number tutrb $CFLT 


enforce zero radial and tangential velocities for post processing 


execute front and back ends only 


employ froude damping 


employ dirt lib over set technology 


employ grid adaption 


exclude continuous field from carver volume fraction nonnalization 


exclude continuous field from carver 


farfield grid smoothing 


farfield patch species mass fractions $VAR 


farfield patch species volume fractions $VAR 


false time step for u $FTU_fieldl, $FTU_fieldl,...,$FTU_fieldN 


false time step for v $FTV_fieldl, $FTV_fieldl,...,$FTV_fieldN 


false time step for w $FTW_fieldl, $FTW_fieldl,...,$FTW_fieldN 


false time step for a $FTA_fieldl, $FTA_fieldl,...,$FTA_fieldN 


false time step for h $FTH_fieldl, $FTH_fieldl,...,$FTH_fieldN 
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false time step for k $FTK field 1, $FTK fieldl, ...,$FTK fieldN 


false time step for ur$FTUU field 1, $FTUU field 1,...,$FTUU fieldN 


false time step for vv $FTVV Field 1, $FTVV field 1,...,$FTVV fieldN 


false time step for ww $FTWW fieldl, $FTWW field 1,...,$FTWW fieldN 


false time step for uv $FTUV_fieldl, $FTUV_fieldl,...,$FTUV_fieldN 


false time step for vw $FTVW fieldl, $FTVW fieldl, ...,$FTVW fieldN 


false time step for wu $FTWU_fieldl, $FTWU_fieldl,...,$FTWU_fieldN 


false time step for f $FTF_fieldl, $FTF_fieldl,...,$FTF_fieldN 


false time step for e $FTE_fieldl, $FTE_fieldl,...,$FTE_fieldN 


fluid rheology $FRHE fieldl, SFRHE field 1,..., SFRHE fieldN 


employ mixture mass for pressure corrector 


employ mixture volume for pressure corrector 


employ cpe continuity coupling 


employ symmetric form of modified pressure 


employ thermodynamic data fits $THFIT_fieldl, $TFIFIT fieldl,. ..,$THFIT_fieldN 


function entry/exit echo on 


function entry/exit echo off 


employ leonard exact diffusion term 


employ gibson modeling 


employ crusader modeliiSCRUELECTRHXEFF, $CRU_TRANSOILHX_EFF, 

$CRU ENGINEHX EFF, $CRU ELECTRHX LOSS, 
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CRUTRAN S OILHX,LO S S , $CRU_ENGINEHX_LOSS 


employ model 5415 modeling 


employ merkle deutsch flate plate modeling 


employ hiplate modeling 


employ darpa 12 inch modeling 


employ large flat plate modeling 


employ sub off modeling 


employ asds modeling 


employ uwx cavitator modeling 


employ mruuv modeling 


employ nrc case modeling 


employ elgho bashi modeling 


employ bubble rise modeling 


employ meghan modeling $IMEGHAN 


employ bistline coding $UCUR, $VCUT, $WCUR, $OXYCUR, $OYZCUR, 

$OZXCUR, 

SXCENTCUR, SYCENTCUR, $ZCENTCUR 


employ haworth lungmodeling 


employ visitor center modeling 


employ cfd ship nphase bubble transport procedure 


gibson modeling starting timestep 
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helicity filter $HELICITY FILTER 


helicity smoothing SHELICITY SMOOTH 


ensight gid not written 


flow velocity initialization using stringer 


hard coded pressure distributions 


hard coded inlet 


hard coded inlet swirl for cyclic verification 


field 1 constant molecular viscosity $VISMC 


field 1 constant density $RHOC 


field 1 constant surface tension $SIGMAC 


free surface model $NFREESURF_surfacel, 

NFREESURF surface2,...,NFREESURF surfaceN 


homogeneous gas mixture fields 


homogeneous gas mixture field parameters $JANNAFOPT 


homogeneous incompressible mixture fields 


homogeneous incompressible mixture field parameters $VAR1, $VAR2 


initialize u 

field $u_fieldO, $u_fieldl,...,$u_fieldN 


initialize v field $v_fieldO, $v_fieldl,...,$v_fieldN 


initialize w 

field $ w_fieldO , $ w_field 1 , . . . , $ w_fieldN 


initialize a 

field $a_fieldO, $a_fieldl,...,$wa_fieldN 
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initialize ai field $ai_fieldO, $ai_fieldl,...,$ai_fieldN 


initialize species mass fractions $ VAR 1, $VAR2 


initialize species volume 


fractions 

$VAR1 


initialize 

k field $K fieldO, $K fieldl,...,$K fieldN 


initialize e field $E fieldO, $E fieldl,...,$E fieldN 


initialize uu field $UU_fieldO, $UU_fieldl,...,$UU_fieldN 


initialize vv field $VV fieldO, $VV fieldl,...,$VV fieldN 


initialize ww 

field $WW fieldO, $WW fieldl,...,$WW fieldN 


initialize uv field $UV_fieldO, $UV_fieldl,...,$UV_fieldN 


initialize vw 

field $VW_fieldO, $VW_fieldl,...,$VW_fieldN 


initialize wu 

field $WU fieldO, $WU fieldl,...,$WU fieldN 


initialize f 

field $F fieldO, $F fieldl,...,$F fieldN 


initialize pressure using s req 


initialize with user initialize field 
routine $VAR1 


inlet grid smoothing 


inlet patch species mass 


fractions 

$MFRAC N 


inlet patch species volume 


fractions 

$VFRAC N 
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interior wall $WPATCH_field(), 

patch SWPATCH fieldl,...,$WPATCH fieldN 


intrinsic swirl 

filter $ISWRL FILTER 


intrinsic swirl 


smoothing 

$ISWRL SMOOTH 


ignore z direction indes delta computation 


initial velocities specified in cylindrical coordinates 


inlet profiles two dimensional specified in inflow.pro 


inlet species profiles specified in inflow.pro 


inlet species profiles 2d specified in inflow.pro 


interfacial area does not feedback into bubble diameter 


interfacial area coalescence model for each field $coales_aint_model_fieldO, 

$coales_aint_model_fieldl , . . . , 
$coales aint model fieldN 


interfield drag models $drag_fielda, $drag_fieldb, $drag_modelid, $drag_usermultiplier 


interfield non drag models $nondrag_fielda, $nondrag_fieldb, $nondrag_modelid, 

$nondrag_usennultiplier, $nondrag_coeff_id2, 


interpolation scheme for face values SFACEINTERP 


incorporate interior wall interface forces $INTERFACE_FORCES 


interpret stringer coordinates as cylindrical 


impose wall bubble diameter kinematic constraint 


immersed boundaries 


kumar bin partitioning $KUMAR_MOMl, $KUMAR_MOM2 


kumar bin characteristic diameter $KUMAR DIA 
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kumar bin characteristic volume $KUMAR VOL 


kumar breakup model $ KUMAR BREAKUP MODELID, $KUMAR BMULT 


kumar coalescence model $KUMAR COAL MODELID, $KUMAR CMULT 


lake evaporation model $EVAP MODELID, $EVAP FACEID 


limit dissipation in production tenn $EPS_PROD_MAX 


limit dissipation in kumar breakup model $EPS_BREAKUP_MAX, 
$LIMIT EPS BREAKUP 


limit er fonn ean flow $LIMIT 


limit er for turbulence $LIMIT2EQ 


inviscid flux scheme $IFLX 


pet sc specify solver iterations $IPETSC_ITER 


pet sc specify solver tolerance $IPETSC_TOL 


pet sc specify solver $IPETSC_SOLVR 


pet sc specify preconditioner SIPETSC PRECON 


pet sc print nonn $IPETSC_PRINTN 


print t min and t max 


modified production kato launder 


moving grid 


moving grid read frequency $GMR_FREQ 


moving grid read grids 


moving grid compute grids 


moving grid compute grids prescribed motion 


multiple frames of reference 
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pad vertex coordinates 


neglect off diagonal terms in pressure poisson equation 


parallel strategy for pressure corrector: matrix level 


parallel strategy for pressure corrector: explicit partition boundary update 


pressures at pressure boundaries specified transiently 


perforin agglomeration operation 


print linear solver residual severys weep 


print linear solver residuals after final sweep 


overwrite density by mid $RHOC_MID 


overwrite molecular viscosity by mid $VISMC_MID 


overwrite gas molecular viscosities $VISGAS_fieldl, $VISGAS_field2,..., 
$VISGAS fieldN 


produce viewable perprocessor ensight output files 


produce ensight files with partition boundaries 


produce ensight files without partition boundaries 


produce ensight files with i blanked cells 


produce ensight files without j blanked cells 


produce ensight output 


produce ensight scalar files for ai 


produce ensight scalar files for a 


produce ensight scalar files for temperature 


produce ensight scalar files for h 
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produce ensigh tscalar files for k 


produce ensight scalar files for e 


produce ensight scalar files for v2f parameters 


produce ensight scalar files for reynolds stresses 


produce ensight scalar files for eddy viscosity 


produce ensight scalar files for density 


produce ensight scalar files for helicity 


produce ensight scalar files for intrinsic swirl 


produce ensight scalar files for debug varO 


produce ensight scalar files for debug varl 


produce ensight scalar files for debug var2 


produce ensight scalar files for debug var3 


produce ensight scalar files for debug var4 


produce ensight scalar files for debug var5 


produce ensight scalar files for debug var6 


produce ensight scalar files for debug var7 


produce ensight scalar files for debug var8 


produce ensight scalar files for debug var9 


produce ensight scalar files for species mass fractions 


produce ensight scalar files for i bla nk 


produce fieldview output 


produce data explorer output 
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produce wall boundary layer output 


pressure gradients req stream sheet correction 


production destruction ratio clip $PRODCLIP 


perfect gas compressibility parameters $RGAS_fieldl, $RGAS_field2,...$RGAS_fieldN 
SSPRATIO field 1, SPRATIO field2,..., SPRATIO fieldN 


calorically perfect gas compressibility parameters $RGAS_fieldl , $RGAS_field2, 

. . . $RGAS_fieldN 

SSPRATIO field 1, SPRATIO field2,..., SPRATIO fieldN 


print pressure boundary reverse flow information 


print outlet boundary reverse flow information 


print cylindrical coordinate patchfiles 


print fieldview in cylindrical coordinates 


print pressure patch profile for inflow profile 


pressure boundary use ambient velocities for inflow 


pressure boundary use extrapolated velocities for inflow 


prandtl number laminar $PRDTLL_fieldl, $PRDTLL_field2,..., SPRDTLL fieldN 


prandtl number turbulent SPRDTLT field 1, SPRDTLT field2, . . . , SPRDTLT fieldN 


mass diffusion coefficient species laminar $LMASSDIFF_fieldl, 
SLMASSDIFF field2,,.., $LMAS SPIFF fieldN 


mass diffusion coefficient species turbulent $TMASSDIFF_fieldl, 
STMASSDIFF field2,..., $TMASSDIFF fieldN 


porous wall patch species mass fractions 


porous wall patch species volume fractions 


pressure grid smoothing 


pressure patch species mass fractions 
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pressure patch species volume fractions 


pressure profile patch 


pid attributes 


number of additional species $NPHI 


order of accuracy of inviscid terms mean flow 


order of accuracy of inviscid terms turbulence 


mass transfer model cprod rprod c destr dest $CPROD, $RPROD, $CDEST, $RDEST 


mass transfer model flag $MTMODEL 


mass transfer model linearization flag $ILMTRANS 


natural cavitatinonumber $CAVNUM 


pre conditioner flag $IPRECON 


pre conditioning parameter beta $BETA 


minimum number of sgs sweeps $ISGSMIN 


maximum number of sgs sweeps $ISGSMAX 


number of sgs sweeps $ISWEEP 


solver choice for velocity components jacobi 


solver choice for velocity components jacobi uvw 


solver choice for velocity components petsc 


solver choice for pressure amggs $SCPRESS_fieldl, $SCPRESS_field2,..., 
$SCPRESS fieldN 


solver choice for pressure amgilu $SCPRESS_fieldl, $SCPRESS_field2,..., 
$SCPRESS fieldN 


use directional coarsening $ISOCOARSE_fieldl, $ISOCOARSE_field2,..., 
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$ISOCOARSE fieldN 


use direct solve on coarsest grid $DIRSOLVE_fieldl, $DIRSOLVE_field2,..., 
$DIRSOLVE fieldN 


solver choice for pressure block correction $NIN, $NJN, $NKN, $PCORR_fieldl, 
SPCORR field2,..., $PCORR fieldN 


solver choice for pressure petsc 


solve mixture momentum mass centered mixture velocity 


solve mixture momentum volume centered mixture velocity 


read wall temperature data from file 


read wall match data from file 


read wall proximity from file 


restart with fieldO frozen 


restart with pressure frozen 


solver monitor iteration $ITER MON 


solver choice for turbulence scalars petsc 


solver choice for all scalars ilu 


solver choice for enthalpy jacobi 


solver choice for enthalpy petsc 


relaxation factor for a $RFA 


relaxation factor for ai $RFAI 


relaxation factor for s $RFS 


relaxation factor for uu $RFUU 


relaxation factor for vv $RFVV 
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relaxation factor for ww $RFWW 


relaxation factor for uv $RFUV 


relaxation factor for vw $RFVW 


relaxation factor for wu $RFWU 


relaxation factor for f $RFF 


relaxation factor for p $RFP 


rhie chow multiplier $RHEICHOW_MULT_fieldl, $RHEICHOW_MULT_field2,..., 
$RHEICHOW MULT fieldN 


surface roughness height $SURF HEIGHT 


residual print file not written 


restart files not written 


transient file write frequency $RESTART_FREQ 


solver sweeps for k SNSWEEP field 1, $NS WEEP fieldl,...,$NSWEEP fieldN 


solver sweeps for e SNSWEEP field 1, $NS WEEP fieldl,...,$NSWEEP fieldN 


solver sweeps for species $NSWEEP_fieldl, $NSWEEP_fieldl,...,$NSWEEP_fieldN 


rectilinear grid 


second order viscous bcs for hexs and prisms 


solve energy equation only 


solve energy equation with frozen flow field 


solve species equation only 


solve species equation with frozen flow field 


relaxation factor applied to linear solver only 


thin layer approximation employed for momentum 
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thin layer approximation not employed for momentum 


thin layer approximation employed for scalars 


thin layer approximation not employed for scalars 


spatial discretization interfacial area density $ISPATIAL 


spatial discretization species equation $ISPATIAL 


spatial discretization density $ISPATIAL 


spatial discretization enthalpy $ISPATIAL 


temporal discretization interfacial area density $ITEMPOR 


temporal discretization species equation $ITEMPOR 


temporal discretization enthalpy $ITEMPOR 


temporal discretization grid $ITEMPOR 


temporal discretization mass $ITEMPOR 


temporal discretization $ITEMPOR 


use absolute velocities as dependent variables 


set pressure correction to zero on partition boundaries 


set pressure correction gradient to zero on partition boundaries 


use block correction to update pressure correction on partitions 


turbulent flow menter k omega 


turbulent flow menter k omega sst 


turbulent flow goldber gkr 


turbulent flow v2f 


turbulent flow frsm 
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turbulent flow wilcox komega 


turbulence model for each field 


turbulence model reynolds number regime for each field 


strongly coupled compressibility 


update enthalpy during pressure corrector 


update species during pressure corrector 


use ficks law form for mass diffusion 


use scalar relaxation for field coupling 


use block relaxation for field coupling 


use scalar relaxation for pc and re coefficients 


use block relaxation for pc and re coefficients 


symmetry grid smoothing 


symmetry patch 


smooth volume fraction in backend $VSMOOTH 


solid region specification mid density conductivity cp $RSQL1D, $KSQL1D, SCSOL1D 


sgs parameters is chkgs tolsor fact $1SCHK, $GSTOL, $SORFACT 


z test, no input yet 


weakly coupled compressibility 


wall grid smoothing 
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Software Delivery and Installation Summary 

The gzipped tar file with which the NPHASE-PSU software is delivered unpacks 
into the following UNIX directory structure: 

tar -xzvf nphase-*****.tar.gz 

setup (bash script to set up UNIX environment for NPHASE-PSU) 

petsc-xxxxx (where xxxxx is the latest version of PETSC delivered with the code, tar 
install file) 

Directories: 

srcbase 

FUMP 

EMERGE 

METIS 

SORT CYCLIC 
SUGGARDIRT 
TUTORIALl 

TUTORIAL HIPLATE (not included in the NGRC distribution) 

TUTORIAL 5415 (not included in the NGRC distribition) 

TUT ORI AL_M ACHB UMP 
TUTORIALGEAR 

Running Tutorials on a System where NPHASE is Already Installed 

In this case delete the all of the directories except the tutorial directories, copy the 
nphase executable, nphase, into each of the tutorial directories, and set the following 
aliases: 

alias fump ’${NPHASE-PSU_PATH} /FUMP/fump' 

alias emergetrans '${NPHASE-PSU_PATH} /EMERGE/emergetrans' 

alias emerge '${NPHASE-PSU_PATH}/EMERGE/emerge' 

alias sort cyclic '${NPHASE-PSU_PATH}/SORT CYCLIC/sort cyclic' 

with the path, ${NPHASE-PSU_PATH} set accordingly. These utilities each need to be 
recompiled as follows: 

1) In directory METIS: tar -zxvf metis-4. 0.3. tar. gz; cd metis-4.0.3; make 

2) In directory FUMP: make 

3) In directory EMERGE: make -f makeem ; make -f makeemt 

4) In directory SORT CYCLIC: make 
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Installing NPHASE-PSU software and Running Tutorials 

The delivered NPHASE-PSU source code resides in the directory srcbase. The code 
unpacks with numerous source fdes (*.c), include fdes (*.h), fdes used for installation 
(makefile. linux, Obj nphase) in that directory. 

To compile NPHASE-PSU, edit the file “setup”, to point the installation to where 
the user wishes NPHASE-PSU to reside on their system. On line 2, modify the 
environment variable $NPHASE_OBJ set to ${NPHASE-PSU_PATH}/build. The same 
modifications need to be carried out for PETSC DIR, DIRTLIBHOME and 
P3DLIB HOME and MPI HOME. The variable PETSC ARCH should be set to reflect 
the compiler architecture used. If using Intel compilers and MPI, the four source 
/software/intel/*** lines should be changed to reflect the location of the users’ Intel 
compiler installation, and the same should be done for LD LIBRARY PATH variables. If 
not using Intel compilers and MPI, the four source statements are unnecessary, and can be 
removed. 

The “setup” script should be called at login whenever NPHASE-PSU is to be used 
to ensure enviromnent variables are set appropriately. 

Compilation of NPHASE-PSU depends upon DiRTlib and P3Dlib. To compile 
these, after making appropriate modification to and calling the setup script, perform the 
following: cd ${P3DLIB_HOME}; ../configure — prefix=$ {P3DLIBHOME}; make; 
make install. Then: cd ${DIRTLIB_HOME}; ../configure — prefix=${DIRTLIB_HOME} - 
-with-p3d=${P3DLIB_HOME} — with-mpi=${MPI_HOME}; make; make install. 
Compilation of DiRTlib with a specified ${MPI_HOME} directory requires that it contain 
${MPI_HOME}/bin, lib, include. If using Intel MPI in 64 bit mode, the directory 
${I_MPI_ROOT} will contain bin64, lib64, include64 subdirectories. To get DiRTlib to 
compile correctly, simply create a dummy directory containing symbolic links bin, lib, 
include that point to these, and point ${MPI_HOME} to this dummy directory within the 
setup script. One caveat to this is that PETSc should be compiled by directing it to the 
actual installation of Intel MPI, as opposed to the dummy one. The procedure for 
accomplishing this is available on the PETSc website. 

Upon appropriate modification to the setup script, and successful installation of 
PETSc, DiRTlib, and P3Dlib, the user can now compile NPHASE-PSU by: 

cd ${NPHASE-PSU_PATH}/srcbase 

make -f makefile.linux -j n 

(where n is the number of processors on the compile node to be used for compilation 
[typically the number of processors on the head node of a cluster]). 

The makefile automatically links to appropriately compiled PETSC libraries 
(PETSC is the open source linear solver library used by NPHASE) and “DIRTLib” 
libraries (DIRTLib is the open source overset meshing library used by NPHASE), compiles 
all of the objects and generates the final executable: nphase. To run the tutorials, copy the 
nphase executable into each of the tutorial directories. 
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